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TECHNICAL  SUMMARY 


The  objective  of  this  project  is  to  determine  the  yield  bias  of  underground  nuclear  tests  induced  by 
the  presence  of  a  high  velocity  descending  slab  beneath  the  test  site.  Specifically,  the  effect  of  the 
Aleutian  slab  is  being  investigated  on  the  US  underground  tests  Longshot,  Milrow,  and  Cannikan. 
P  wave  seismograms  will  be  synthesized  using  dynamic  ray  tracing  and  superposition  of  Gaussian 
beams  in  three-dimensional  models  of  the  Aleutian  slab  determined  from  P  travel  time  delays. 
Focusing  and  defocusing  and  multipathing  at  teleseismic  distances  will  be  evaluated  by  comparison 
of  observed  with  synthetic  seismograms  of  the  Aleutian  tests. 

To  calculate  effects  of  focusing  and  defocusing  of  three  dimensional  slab  structure  an  alternative 
method  of  dynamic  ray  tracing  and  Gaussian  beam  superposition  was  developed,  requiring  a  fewer 
number  of  equations  to  be  integrated  along  ray  and  eliminating  the  need  to  evaluate  second  spatial 
derivatives  of  velocity.  We  term  this  method  vicinity  ray  tracing.  Its  derivation  and  application  is 
described  in  a  reprint,  which  is  included  as  the  section  I  of  this  report. 

This  method  was  applied  to  predict  the  amplitudes  of  P  waves  radiated  by  nuclear  tests  in  the 
Aleutians  using  three  dimensional  models  of  the  Aleutian  slab  proposed  from  the  fitting  of  P  travel 
time  residuals  of  Aleutian  seismic  events.  Section  II  compares  these  synthetic  amplitudes  with 
observed  distributions  of  m*,  residuals  of  Aleutian  tests.  The  results  demonstrate  that  a  network 
average  of  mi  measured  from  stations  concentrated  in  the  shadow  zone  of  the  Aleutian  slab  can 
lead  to  an  underestimate  of  the  size  of  seismic  events  in  the  Aleutian  island  ridge  by  up  to  OAnii 
units  compared  to  the  same  size  events  located  in  regions  unaffected  by  focusing  and  defocusing. 
This  underestimate  is  reduced  to  0.1  mi  units  for  the  most  probable  distribution  of  stations  in  a 
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network  average  of  m*.  Azimuthal  weighting  of  stations  in  such  an  average  will  not  significantly 
reduce  this  underestimate  because  of  the  large  areal  extent  of  the  slab  shadow  zone  and  its  overlap 
with  a  large  concentration  of  European  and  Canadian  seismic  stations.  The  need  for  a  relative 
correction  for  focusing  and  defocusing  effects  in  comparisons  of  NTS  and  Aleutian  tests,  however, 
may  net  be  justified  due  to  the  documented  presence  of  a  defocusing  structure  beneath  at  least  a 
portion  NTS  that  can  produce  a  similar  sized  negative  mi  bias. 

Section  III  is  preprint  of  a  paper  accepted  for  publication,  which  treats  a  problem  important 
to  nuclear  monitoring  at  regional  distance.  The  locked  mode  method  of  synthesizing  complete 
regional  seismograms  was  modified  to  include  the  Langer  uniform  asymptotic  approximation  to 
vertical  wavefunctions  within  layers  having  linear  vertical  velocity  gradients.  Computational  ex¬ 
periments  were  made  to  (1)  quantify  the  breakdown  in  the  asymptotic  approximation  to  the  vertical 
wavefunctions  as  frequency  decreases  am  /or  magnitude  of  the  vertical  gradient  increases  and  to 
(2)  illustrate  and  review  some  of  the  elfects  of  the  vertical  velocity  gradients  on  the  propagation  of 
regional  phases.  This  preprint  is  an  update  of  the  paper  included  in  the  third  technical  report  of 
this  project.  It  includes  additional  examples  and  tests  not  shown  in  that  report  as  well  corrections 
of  equations  given  in  appendices.  ltd  conclusions  are  unchanged  from  those  given  in  the  earlier 
draft. 
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SUMMARY 

A  method  for  computing  seismic  wavefieids  in  a  high-frequency  approximation  is 
proposed  based  on  the  integration  of  the  kinematic  ray  tracing  equations  and  a  new 
set  of  differential  equations  for  the  dynamic  properties  of  the  wavefront,  which  we 
call  the  vicinity  ray  tracing  (VRT)  equations.  These  equations  are  directly  obtained 
from  the  Hamiltonian  in  ray  centred  coordinates,  using  no  paraxial  approximations. 
This  system  is  comparable  to  the  standard  dynamic  ray  tracing  (DRT)  system,  but  it 
is  specified  by  fewer  equations  (four  versus  eight  in  3-D)  and  only  requires  the 
specification  of  velocity  and  its  first  spatial  derivative  along  a  ray.  The  VRT 
equations  describe  the  trajectory  of  a  ray  in  the  ray  centred  coordinates  of  a 
reference  ray.  Quantities  obtained  from  vicinity  ray  tracing  can  be  used  to 
determine  wavefront  curvature,  geometric  spreading,  traveltime  to  a  receiver  near 
the  reference  ray,  and  the  KMAH  index  of  the  reference  ray  with  greater  numerical 
precision  than  is  possible  by  differencing  kinematically  traced  rays.  Since  second 
spatial  derivatives  of  velocity  are  not  required  by  the  new  technique,  parametriza- 
tion  of  the  medium  is  simplified,  and  reflection  and  transmission  of  beams  can  be 
calculated  by  applying  Snell’s  law  to  both  vicinity  rays  and  central  rays.  Conversion 
relations  between  VRT  and  DRT  can  be  used  to  determine  the  paraxial  vicinity  of 
DRT,  in  which  the  errors  of  the  paraxial  approximations  of  DRT  remain  small.  In 
cither  DRT  or  VRT,  the  width  of  Gaussian  beams  can  be  physically  defined  from 
the  width  of  the  Fresnel  volume  surrounding  the  central  ray.  Because  no  paraxial 
approximations  are  made,  the  superposition  of  the  Gaussian  beams  defined  from 
vicinity  rays  should  exhibit  a  much  slower  breakdown  in  accuracy  as  the  scale  length 
of  the  medium  given  by  u/|Vn|  approaches  the  beamwidth. 

Key  words:  asymptotic  ray  theory,  dynamic  ray  tracing,  seismic  wavefieids. 


i  INTRODUCTION 

Many  high-frequency  asymptotic  solutions  of  the  wave 
equation  have  been  developed  as  effective  tools  for 
computing  wave  fields  in  inhomogeneous  3-D  media.  Two  of 
the  most  widely  applied  are  the  WKBJ/Maslov  method 
(Chapman  1978.  Chapman  &  Drummond  1982)  and  the 
Gaussian  beam  method  (Babii  &  Kirpiimkova  1979,  Popov 
1982;  Cerveny,  Popov  &  PSeniik  1982;  Cerveny  &  PSeniik 
1983)  Both  of  these  techniques  cstur.alt  the  kinematic  and 
dynamic  properties  of  a  wavefront  from  approximate 
solution;,  to  the  elastodynamic  wave  eq^atiun  based  on  ray 
theory  The  superposition  techniques  of  Gaussian  beams 
and  WKBJ  plane  w..ves  as  well  as  their  stationary  phase 
approximation  m  geometric  tay  theory  ail  require  similar 
amplitude  and  weighting  functions.  These  amplitude 
function'  can  l*.  found  by  integrating  a  system  of  equations 
known  as  the  dynamic  ray  tracing  (DRT)  equations. 


The  DRT  equations  can  be  derived  from  cither  the 
eikonal  equation  by  substitution  of  a  paraxial  approximation 
(Cerveny  &  Hron  1980;  Cerveny  1985),  or  from  the 
parabolic  wave  equation  (Cerveny  el  al.  1982;  Popov  1982; 
Cerveny  &  PSenftk  1983).  The  DRT  equations  have  both 
limitations  and  complications.  The  limitations  are  associated 
with  the  use  of  the  paraxial  approximation,  and  the 
complications  are  due  to  the  use  of  multiple  coordinate 
systems. 

Tbs  limitations  associated  with  the  paraxial  approxima¬ 
tion  are  exhibited  whenever  the  DRT  equations  are  used  to 
estimate  the  traveltime  and  amplitude  at  a  point  in  the 
neighbourhood  of  a  ray  from  a  second-order  Taylor 
expansion  of  the  wavefront  at  a  point  along  the  ray.  The 
Tayloi  expansion  is  the  essential  step  in  the  definition  of 
Gaussian  beams  an  paraxial  rays.  The  region  in  which  the 
error  of  this  Taylor  expansion  remains  below  some  specified 
threshold  is  generally  referred  to  as  the  paraxial  vicinity. 
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The  fundamental  problem  with  the  paraxial  approximation 
is  that  it  is  not  simple  to  specify  the  spatial  bounds  of  the 
paraxial  vicinity  in  a  three-dimensionally  varying  model.  In 
general,  one  must  not  attempt  to  evaluate  the  Taylor 
expansion  too  far  from  the  central  ray.  but  it  is  unknown 
how  the  error  grows  away  from  the  central  ray. 

The  complications  associated  with  the  use  of  two 
coordinate  systems  are  best  appreciated  by  considering  the 
most  general  case  of  a  three-dimensionally  varying  medium. 
In  three-dimensionally  varying  media,  the  usual  approach  is 
to  specify  the  DRT  equations  using  two  coordinate  systems: 
ray  coordinates,  usually  consisting  of  the  take-off 
angle/azimuth  at  the  source,  and  ray  centred  coordinates, 
consisting  of  an  orthogonal  curvilinear  system  that  moves 
along  with  the  ray  (Fig.  1),  The  use  of  two  coordinate 
systems,  while  having  the  advantage  of  converting  a 
non-linear  Ricatti  equation  into  a  system  of  linear 
equations,  increases  the  number  of  equations  needed  to 
describe  the  quantities  affecting  the  amplitude  of  the 
wavefield.  In  either  the  fixed  Cartesian  or  ray  centred 
coordinates,  the  standard  DRT  equations  require  the 
specification  cf  the  second  spatial  derivatives  of  velocity 
along  a  ray.  This  either  forces  the  model  to  be  parametrized 
with  continuous  first  derivatives  of  velocity  or  complicates 
the  integration  by  requiring  jump  conditions  on  the  dynamic 
quantities.  These  jump  conditions  are  obtained  by  matching 
the  paraxially  estimated  phase  on  either  side  of  a 
discontinuity  in  gradient  (Ccrveny  &  Hron  1980;  Ccrveny 
1985). 

In  this  paper,  we  develop  alternative  methods  for 
calculating  wavefront  curvature,  geometric  spreading,  and 
the  widths  of  Gaussian  beams  that  eliminate  many  of  these 
problems.  These  alternative  methods  are  based  on 
quantities  calculated  from  a  system  of  equations  for  the  path 
of  a  ray  near  a  reference  ray  (Fig.  2).  This  nearby  or 
‘vicinity’  ray  may  be  calculated  exactly  by  the  kinematic  ray 
tracing  equations,  with  geometric  spreading  and  wavefront 
curvature  estimated  by  ray  differencing  (Gajewski  & 
PSentik  1987),  or  determined  approximately  by  integrating  a 
system  of  equations  we  refer  to  as  the  'vicinity  ray  tracing 
system  tVRT)’.  Gaussian  beams  arc  defined  by  assigning  a 
Gaussian  distribution  of  amplitude  to  each  central  ray.  The 
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Figure  1.  The  ray  centred  coordinates  (s,  I  is  the  unit 

tangent  vector  of  a  central  ray  and  e,  and  e,  are  the  unit  normal 
vectors  to  l.  The  coordinate  s  measures  the  arclength  along  a 
central  ray  from  an  arbitrary  reference  point.  qt  and  q2  represent 
length  coordinates  and  form  a  2-D  orthogonal  coordinate  system  in 
the  plane  normal  to  J2  at  O. 


C| 


Figure  2.  The  geometry  of  the  vicinity  ray  tracing  system:  rj,  is  the 
angle  difference  between  the  tangential  vector  ?  of  a  central  ray  and 
the  tangential  vector  of  a  vicinity  ray  in  the  i-e,  plane,  q,  is  the 
distance  between  the  central  ray  and  the  vicinity  ray  in  the  t-e, 
plane  at  S. 

beamwidth  of  the  Gaussian  is  taken  to  be  the  half-width  of 
the  beam  Fresnel  volume  surrounding  the  central  ray.  The 
outermost  vicinity  ray  in  a  beam  Fresnel  volume  has  a 
half-period  traveltime  difference  with  respect  to  that  of  the 
central  ray.  Since  beamwidths  are  related  to  the  beam 
Fresnel  volume,  diffracted  wavefields  can  be  accurately 
estimated  by  a  superposition  of  Gaussian  beams  without  the 
ambiguity  associated  with  a  freely  varying  beamwidth 
parameter.  The  geometrical  spreading  of  a  wavefront  is 
computed  from  the  conservation  law  of  energy  in  VRT, 
rather  than  from  the  direct  solution  of  the  transport 
equation  as  in  DRT. 

In  the  following  sections  we  first  review  the  derivation  of 
the  standard  DRT  system  and  the  limitations  of  the  paraxial 
approximation.  Next  the  vicinity  ray  tracing  system  is 
derived  from  the  Hamiltonian,  in  which  no  paraxial 
approximations  are  made.  Expressions  for  the  traveltime 
and  wavefront  curvature  in  the  neighbourhood  of  a  central 
ray  are  derived  using  this  system.  Gaussian  beams  arc 
defined  using  vicinity  rays  to  approximate  the  beam  Fresnel 
volume.  The  accuracy  of  estimated  geometrical  spreading, 
ray  trajectory,  and  the  traveltimes  of  vicinity  rays  is  tested 
and  compared  in  a  simple  1-D  model  using  VRT  and  DRT. 
Finally,  seismograms  are  synthesized  and  compared  in  the 
same  model  using  a  superposition  of  Gaussian  beams 
derived  from  VRT  and  the  WKBJ  method. 

2  A  REVIEW  OF  DYNAMIC  RAY  TRACING 
(DRT) 

2.1  Ray  centred  coordinates 

Consider  an  arbitrary  ray  corresponding  to  a  P- wave  ahd 
introduce  ray  centred  coordinates  s,  qx,  q2  (Fig-  0-  The 
orthogonal  ray  centred  coordinate  system  along  the  central 
ray  Q  and  its  computations  are  described  in  Popov  & 
PSenCik  (1978),  and  Cerveny  &  Hron  (1980).  The  ray 
centred  coordinates  are  limited  to  a  vicinity  of  the  origin 
(q,  =0)  in  which  the  central  ray  field  is  regular.  In  Fig  1, 
the  coordinate  j  measures  the  arclength  along  a  central  ray 
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from  an  arbitrary  reference  point.  <7,  and  q2  represent 
length  coordinates  and  form  a  2-D  Cartesian  coordinate 
system  in  the  plane  normal  to  Q  at  O,  with  origin  at  £2.  All 
three  components  (s,  q„  q2)  in  the  ray  centred  coordinate 
system  depend  on  the  azimuth  and  vertical  take-off  angle 
( <p ,  6).  The  basis  of  the  coordinate  system  forms  a 
right-handed  system  of  the  three  unit  vectors  i,  e,  and  e2 
where  i  is  the  unit  tangent  vector  to  the  central  ray  Q. 


2.2  The  paraxial  approximation  and  its  limitations 

The  standard  DRT  system  can  be  derived  from  either  the 
eikonal  equation  (Cerveny  &  Hron  1980;  Madariaga  1984; 
Cerveny  1985)  or  from  the  parabolic  wave  equation  (Popov 
1982;  Cerveny  &  PSenCik  1983).  In  either  derivation,  a 
paraxial  approximation  is  assumed  at  some  stage,  which 
involves  a  Taylor  expansion  of  the  waveficld  about  the 
central  ray.  The  terms  are  omitted  in  approximations  often 
without  specifying  validity  conditions. 

To  illustrate  the  problems  with  the  DRT  system,  let  us 
review  the  derivation  of  the  2-D  DRT  equations  starting 
from  the  eikonal  equation.  The  eikonal  equation  in  2-D  is 


A  (ilY  (il\7 

h2\ds)  +\dq) 


where  V  =  v(s,  q).  h  is  a  scale  factor  in  the  ray  centred 
coordinates  and  will  be  discussed  subsequently.  The 
traveltime  of  a  vicinity  ray  r (s,  q)  can  be  approximated  at 
q  =  0  (Cerveny  &  Hron  1980;  Cerveny  &  PScnCik  1983; 
Cerveny  1985)  by: 

r(s,  q)  =  r(s)  +  {M{s)q2  (2) 

where  dr/dq  =  0  and  M  =  d2x (s,  q)ldq2.  From  equation 
(2),  it  follows  that 


dx (s,  q) 
dq 


=  Mq, 


(3) 


where  v  =  v(s,  0).  Substituting  equation  (3)  into  (1)  and 
neglecting  higher  order  terms  gives 


M,  +  M2)q2  =  -J2 


1 

v2h2' 


(4) 


The  right  side  of  equation  (4)  can  be  approximated  by 
expanding  the  velocity  V'  up  to  second-order  terms  with 
respect  to  v, 


_1_  1  1 
V2~ v2h2=a  V-qq<l 


(5) 


(Cerveny  &  Hron  1980;  Cerveny  &  PSeniik  1983;  Cerveny 
1985),  where  v  qq  =  d1 vldq 2.  The  standard  DRT  system  is 
obtained  from  equation  (4)  by  using  equation  (5)  and 
expanding  up  to  second  order.  This  gives 

^-+vM2  +  ^-f  =  0.  (6) 

as  v 

Since  the  derivation  of  the  DRT  system  includes 
second  ordei  terms,  any  omitted  terms  must  be  carefully 
evaluated.  Consider  the  scale  factor  h.  The  scale  factor  h  is 


Vicinity  ray  tracing 
given  by  (Cerveny  &  Hron  1980) 

h  =  l+^q.  (7) 

v  v  ' 

Because  the  first-order  term  of  h  is  neglected  in  equation 
(4),  the  term  (2 v  q/v)q  of  /r  in  equation  (4)  must  be 
vanishly  small,  i.e., 

(8) 

The  condition  in  equation  (8)  describes  the  applicability  of 
the  DRT  system.  It  says  that  extrapolation  of  the  wavefield 
at  a  distance  q,  away  from  a  central  ray  using  the  paraxial 
approximation  will  break  down  rapidly  as  the  scale  length  of 
the  medium  decreases,  where  the  scale  length  /  is  defined  by 
/  =  u/|Vu|.  The  extrapolation  distance  must  be  much  less 
than  the  scale  length  of  the  medium.  For  Gaussian  beams,  it 
implies  that  the  beam  width  must  be  much  less  than  the 
scale  length  of  the  medium.  This  can  be  a  severe  restriction 
in  rapidly  varying  models,  in  which  the  criterion  for  validity 
of  ray  theory  (wavelength  «  scale  length)  is  still  well 
satisfied. 


2.3  The  P  and  Q  matrices 

Equation  (6)  is  a  non-linear  ordinary  differential  equation  of 
the  first-order  Ricatti  type.  This  equation  can  be  solved  by 
elementary  analytical  methods.  Following  Cerveny  &  Hron 
(1980),  the  2-D  system  given  by  equation  (6)  can  be 
generalized  to  a  3-D  system  of  linear  differential  equations 
by  introducing  a  2  X  2  matrix  M: 

M  =  U"^Q~1  (9) 
where  Q  is  a  2  x  2  matrix.  Define  a  2  x  2  matrix  P  as 


P  =  iT' 


da 

ds  ' 


(10) 


Substituting  equations  (9)  and  (10)  into  equation  (4),  the 
dynamic  ray  tracing  equations  in  3-D  can  be  written  as 


da  „  dP  1 

—  =uP,  —  =  — iSQ, 
ds  ds  v 


01) 


where  Qq  =  dqjdy,,  P,,  =  dpjdy,  and  y,  are  ray  parameters 
(usually  take-off  angles).  S  is  given  as 


The  DRT  system  has  eight  equations  for  real  Q  and  P, 
arid  16  for  complex  Q  and  P  in  3-D  and  is  specified  in  ray 
centred  coordinates  (s,  qx,  q2)  and  ray  coordinates  (y, ,  y2). 
Cerveny  (1985)  has  shown  that  only  eight  equations  are 
generally  needed  for  Gaussin  beams.  The  DRT  system 
generally  will  have  off-diagonal  terms  in  the  matrices  Q  and 
P.  The  existence  of  these  off-diagonal  terms  is  due  to  the  use 
of  two  coordinate  systems  in  describing  the  equations.  The 
number  of  equations  can  be  reduced  further  if  only  one 
coordinate  system  could  be  used. 
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2.4  Gaussian  beams 

Approximate  and  non-singuiar  solutions  of  the  elastic 
wavelield  can  be  found  in  the  vicinity  of  caustics  and  regions 
of  diffraction  by  superposing  Gaussian  beams  constructed 
from  complex  Q  and  P  matrices.  A  proper  mathematical  or 
physical  meaning  of  complex  parameters  in  the  Q  and  P 
matrices,  however,  is  not  usually  considered  in  routine 
applications  of  the  method.  Complex  Q  and  P  can  be  shown 
to  be  a  consequence  of  an  approximate  solution  for  complex 
rays  emanating  from  a  source  having  a  small  imaginary  part 
to  its  location  in  space  (Felscn  1984;  Wu  1985).  In  practice, 
beamwidths  are  defined  somewhat  arbitrarily  and  are 
adjusted  to  minimize  errors  in  the  beam  superposition 
(KlimcS  1988;  Kim  &  Garmany  1985)  or  tuned  to  minimize 
errors  associated  with  rapid  variations  in  velocity  (Weber 
1988).  White  et  al.  (1987)  have  shown  that  optimum 
beamwidths  strongly  depend  on  the  specific  wave  propaga¬ 
tion  problem  and  the  particular  type  of  boundary 
interactions  occurring  in  the  problem.  One  of  the  reasons 
why  the  concept  of  optimum  beamwidths  does  not  work  well 
is  that  the  energy  of  each  beam  differs  for  different  initial 
beamwidths.  This  is  true  for  all  of  the  various  optimal 
beamwidths  that  have  been  proposed.  If  energy  flux  is  to  be 
conserved  within  a  ray  tube,  then  a  normalization  condition 
must  be  applied  with  respect  to  the  different  initial 
beamwidths.  The  following  section,  which  introduces  the 
vicinity  ray  tracing  system,  discusses  how  such  a 
normalization  condition  can  be  implemented  and  how  beam 
width  in  both  the  DRT  and  VRT  systems  can  be  physically 
related  the  Fresnel  zone  surrounding  a  central  ray. 

3  THE  VICINITY  RAY  TRACING  (VRT) 
SYSTEM 

3.1  Derivation 

Let  us  consider  the  high-frequency  asymptotic  solution  to 
the  wave  equation  in  an  inhomogeneous  medium.  In  order 
to  obtain  the  desired  approximation,  let  us  assume  that  the 
Fourier  component  of  displacement  u  for  frequency  w  is 
expressed  in  the  following  form  in  a  generalized  coordinate 
system: 

u(q„  a)  =  A(q,)c  (12) 

where  i  =  1,2,.  n  is  an  n -dimensional  configuration  space 
whose  coordinates  are  the  n  generalized  coordinates  q,.  A 
and  r  are  an  amplitude  function  and  a  phase  function 
respectively.  Both  A  and  r  are  functions  which  can  be 
assumed  to  be  siowly  varying  with  respect  to  the  wavelength 
A.  The  Hamiltonian  and  the  conservation  law  of  energy  flux 
along  the  wavefron'  aie  applied  in  this  study  to  determine 
the  functions  A  and  r. 

Since  Fermat's  principle  of  least  time  can  be  expressed  by 
the  opciation.1  of  variational  calculus,  a  Lagrangian  and  the 
Hamiltonian  can  be  defined  similar  to  those  used  in 
describing  the  mechanics  of  particles.  For  example,  Farra  & 
Madariaga  (1987)  used  the  Hamiltonian  to  develop  a 
perturbation  theory  to  compute  the  amplitude  and 
traveltime  of  a  vicinity  ray  with  respect  to  a  reference  ray. 

By  applying  Fermat  s  principle,  the  traveltime  t  from 
source  ls,,j  to  receiver  |s,)  can  be  written  as  a  path  integral 


over  the  Langrangian,  L,  by 

*=  f  L(q„qi)ds,  (13) 

•Vo 

where  q^dqjds.  The  Hamiltonian  in  ray  centred 
coordinates  can  be  obtained  from  the  Lagrangian.  The 
Lagrangian  in  ray  centred  coordinates  is  given  as 

i ,  <?2>  <7i.  9a.  s)  =  £ (h2  +  q\  +  q\)xri,  (14) 

where 

h=hx  =  \  +  YJ~ql  (15) 

i«l  v 

and  where  q,  =  dqi/ds,  v.  —  dv/dqi  and  h2  =  h3=l. 
V  =  t/(s,  qt,  q2)  is  the  velocity  of  a  vicinity  ray  and 
u  =  u(s,  0,  0)  is  the  velocity  of  a  central  ray.  The  Lagrangian 
in  equation  (14)  has  s  as  an  independent  variable.  The 
conjugate  momentum  p,  can  be  expressed  as 

p'=W\=vVi F+tr+ti’ 

(16) 


dq2  Vyh2  +  q2  +  ql' 


Equations  (16)  can  be  solved  for  q,  and  q2,  yielding 


Vhpi 


The  Hamiltonian,  H,  is  expressed  as  follows: 

tf(<7,,  qi<  P\>  Pi>s)~  P\q\  +P29z 

-  Uqx,q2,quqi<s)-  (18) 


By  substituting  equations  (17)  into  equation  (18),  the 
Hamiltonian  in  ray  centered  coordinates  is  obtained; 

H(quq2.Pup2.s)=  ~p[l  -  V2(p]+pl)]u2.  (19) 


The  cikonal  equation  can  be  derived  from  the 
Hamiltonian  by  using  the  Hamilton-Jacoby  partial 
differential  equation  with  respect  to  the  arclength  s 
(Goldsten  1980): 


dx 

ds 


q  i>  q  2 


dr  dx 
'  9q dq2 


(20) 


Squaring  both  sides  of  equation  (20)  produces  the  cikonal 
equation  in  ray  centred  coordinates: 


(21) 


The  VRT  system  in  ray  centred  coordinates  can  be 
described  in  terms  of  the  canonical  equations  from  the 
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Hamiltonian  defined  in  equation  (191: 

da i  _  3H  _  hVp, 
ds  3p,  E  ’ 

dq,  J H  _hVp, 
ds  3pz  E 

±1-  LJl 

V  V2E  ’ 


*22) 


ds 


-9 1 


4Pz  =  SamhX£  V,h 

ds  3a-  V  l'*£ " 


where 


E  =  \n-V-(P-~px). 


*23) 


Using  the  cikonal  equation  (21)  in  ray  centred 
coordinates,  the  quantity  £  can  be  approximated  as 


V 

E~'hv 

Substituting  (24)  into  (23)  gives 

dq,  dq,  . 

-  =  vhPi.  —  =  vhP:. 

dp,  h ...  vV  ...hi1 


(24) 


1  —  ^  . 
ds  hv 


V' 


(25) 


dps  h„  vVa.h- 
ds  hv  V'  ' 


Equations  (25)  described  in  terms  of  q,  and  p,  are 
comparable  to  the  dynamic  ray  tracing  equations  (11).  but 
no  paraxial  approximations  have  been  made.  Cerveny  & 
P5cncik  (1979)  derived  similar  equations  to  equations  (25) 
from  the  cikonal  equation.  Cerveny  (1987)  has  also  briefly 
described  the  derivation  of  equations  (25).  showing  how 
they  become  the  standard  DRT  equations  if  paraxial 
approximations  are  substituted  for  h  and  V. 


3.2  Implementation  of  VRT  and  its  relations  to 
kinematic  ray  tracing  and  DRT 

3. 2. 1  The  standard  form 

Cerveny  (1987)  suggested  that  the  system  of  equations  (25) 
would  not  be  suitable  for  calculations  because  the  right-hand 
sides  of  the  equations  for  p,  are  the  difference  of  two  large, 
nearly  equal  quantities.  Note  for  the  case  of  a  vicinity  ray 
very  close  to  the  central  ray  and/or  for  the  special  case  of  a 
medium  that  has  only  very  weak  inhomogencity  the  two 
terms  on  the  right-hand  sides  of  equations  (25)  both 
approach  a  value  equal  to  v  q/v2  In  practice,  we  have 
found  this  not  to  be  a  problem  for  the  test  structure  shown 
in  this  paper  as  well  as  for  most  calculations  involving 
Gaussian  beams  in  media  having  spatial  gradients.  The 
extent  to  which  a  numerical  problem  exists  will  depend  on 
the  structure  and  ray  geometry,  '*  system  (25)  is  used  very 
close  to  a  central  ray  and  the  sp—ial  gradients  are  small, 
then  a  numerical  problem  can  exist.  If  the  system  (25)  is 
used  in  media  having  strong  gradients  and  at  distances 
sufficiently  far  from  the  central  ray.  there  will  not  be  a 


ccss craa3  pio&lsa.  These  are  exarrty  the  skmduses  where 
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3.2.2  A  numerically  ssiblc  form 

An  altercarive  form  of  cqtntxxx  (25)  ezn  be  tsed  in  which  a 
srxxxb  traashioa  can  be  made  besweea  a  form  t irzi  is 
acainte  at  targe  drsraaee s  a  from  the  cmira)  ray  aad  a  fora 
that  remains  carnericalTr  regular  a:  small  dsaaces  from  the 
central  ray  and/or  in  media  having:  very  weak  ta- 
honsogeneity.  This  alternative  form  is  derived  by  describing 
equations  (25)  in  terms  of  q.  2nd  the  aagabr  difference 
between  the  central  ray  aad  the  vicinity  rays.  Let  cs  define 
7j,  as  the  angle  difference  between  the  tangential  vectors  of 
a  central  ray  and  a  vsrimty  ray  in  the  r-e,  phac  and  17.  as 
the  angle  difference  between  the  tangential  vectors  erf  a 
central  ray  and  of  a  vicinity  ray  in  the  i-c,  p!aae  in  ray 
centred  coordinates  (Fig.  2).  q,  is  the  distance  from  the 
central  ray  to  the  vicinity  ray  along  the  e,. 

Let  us  consider  two  specified  vicinity  rays  which  are 
located  in  the  r-e,  plane.  Because  the  r-e,  and  the  r-e. 
planes  are  orthogonal  to  each  other,  q.  and  p.  are  zero  in 
the  former  plane  and  q,  and  p,  are  zero  in  the  latter  plane. 
The  slowness  p,  in  the  r-e,  plane  is  described  as 


_  sin  17, 
Pi  y  • 
r* 


(26) 


where 

K,  =  v(s.  q„  0)  and  V2  =  v(s.  0.  <7.). 

By  substituting  equation  (26)  into  equations  (25).  dqjds  can 
be  rewritten  as 


hfu  . 

-^-=— -sttwj,. 
dr  K 


(27) 


Differentiating  equation  (26)  with  respect  to  s  in  the  ray 
centred  coordinates  yields 


dp,  cos  tj/di],  h/V,^  sin  17, 
ds  V,  ds  V? 


(28) 


where 

V^dVJh.ds. 

dqjds  and  dqjds  can  be  obtained  by  equating  equations 
(25)  and  (28).  Equations  (28)  can  then  be  rewritten  in  the 
following  form  in  terms  of  q,  and  q,  in  ray  centred 
coordinates: 


dq,  h\v  . 

—  =  —  s\nq„ 
ds  tj 


dfh 

ds 


hi v  . 

-7  = -7-sm  77,, 


dq, 

ds 


h,V ,, 


iant),  +  C 


V, 

COS  77,  ' 


dih 

ds 


h,V2, 

V, 


tan  t72  +  D 


Vs 

COS  IJj  ' 


(29) 
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where 


C 


2!S**I 

v]  ’ 


D 


vl  ■ 


traveltime  to  an  arbitrary  receiver  point  can  be  estimated 
from  the  values  of  q t  and  rj,  of  vicinity  rays  near  the 
receiver.  Section  3.5  describes  how  geometric  spreading  can 
be  calculated  using  VRT  without  a  transformation  between 
ray  centred  and  ray  coordinates  by  using  the  conservation 
law  of  energy  flux  along  the  wavefront.  This  section  also 
describes  superposition  of  Gaussian  beams  from  using  VRT. 


If  q  is  very  small  and/or  heterogeneity  is  very  weak,  the 
two  terms  in  the  expressions  for  C  and  D  nearly  equal  one 
another.  In  these  situations,  the  first  term  in  equations  (29) 
is  much  larger  than  the  terms  containing  C  and  D.  Thus,  any 
potential  problem  in  numerical  precision  can  be  easily 
controlled  by  setting  C  and  D  to  zero  in  equations  (29) 
whenever  the  two  terms  defining  C  and  D  equal  one  another 
to  within  three  or  more  significant  digits.  Test  calculations  in 
structures  having  spatial  gradients  with  a  spacing  of  10  or 
more  vicinity  rays  per  medium  scale  length  always  gave  a  C 
and  D  sufficiently  large  that  no  such  check  for  loss  of 
precision  was  necessary. 


3. 2. 3  Dynamic  properties  from  VRT 

From  Fig.  3,  it  is  seen  that  the  curvature  ( K ,)  of  the 
wavefront  at  the  point  (r,  q,)  in  the  t-e,  plane  is  a  function 
of  tan  i},  and  q,\ 


K,= 


tan  1 1, 
9. 


1 


(30) 


Because  the  radius  of  curvature  of  the  wavefront  is 
described  in  terms  of  q,  and  r\,  in  the  VRT  system,  writing 
the  VRT  system  in  terms  of  q,  and  q,  describes  the 
wavefront  of  the  ray  more  directly  thai.  ’he  standard  VRT 
system  written  in  terms  of  q,  and  pr 
Given  the  expression  equation  (30)  for  the  radius  of 
curvature  of  a  wavefront  as  a  function  of  q,  and  q,,  the  VRT 
system  can  be  used  to  estimate  the  traveltime  to  a  point  near 
a  central  ray,  iteratively  solve  the  two  point  ray  tracing 
problem,  and  calculate  geometric  spreading.  To  determine 
all  of  the  relevant  properties  of  the  wavefield,  one  must 
integrate  both  the  kinematic  ray  tracing  system  and  the 
vicinity  ray  tracing  system.  Section  3.3  describes  how 
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Figure  J.  Ray  centred  coordinates  ($.  q,)  and  the  wavefront 
coordinates  (j,  q ,):  the  Jacobian  J  between  two  coordinates  is  given 
in  equation  (45)  The  curvature  K,  and  the  radius  of  the  curvature 
of  the  wavefront  R,  in  the  l-e,  plane  are  described  q[  is  the  normal 
distance  between  0  and  0". 


3.2.4  Dynamic  properties  from  kinematic  ray  tracing 

It  is  worth  noting  that  it  is  possible  to  estimate  most  of  the 
quantities  (29)  from  the  kinematic  ray  tracing  equations. 
The  VRT  system  essentially  traces  a  ray  in  the  ray  centred 
coordinate  system  of  a  reference  ray.  Thus,  in  some 
situations  it  may  be  both  simpler  and  more  accurate  to  trace 
a  nearby  ray  by  integrating  the  kinematic  ray  tracing 
equations  with  a  perturbation  in  the  take-off  angles  of  the 
reference  ray.  The  quantities  <7,  and  q2  can  be  estimated 
from  the  perpendicular  distances  between  one  reference  ray 
and  two  nearby  rays.  Likewise  the  angles  q ,  and  q2  can  be 
estimated  from  the  difference  in  the  local  tangent  of  one 
reference  ray  and  the  local  tangent  of  two  nearby  rays.  The 
loci  of  all  three  of  these  rays  can  be  found  by  integrating  the 
simple  kinematic  ray  tracing  equations  rather  than  the 
vicinity  ray  tracing  equations  given  in  equations  (29). 

In  a  three  dimensionally  varying  medium,  the  two  nearby 
rays  will  not  generally  lie  along  the  directions  of  the  vector 
basis  of  a  ray  centred  coordinate  system  (e,,  e2),  even  if  the 
initial  conditions  were  carefully  chosen  to  achieve  this  result 
at  the  starting  point  of  integration.  Thus,  some  additional 
work  would  usually  be  necessary  to  convert  these 
measurements  into  distances  q,  and  angles  q,  of  rays  lying 
along  the  e,  and  e2  directions. 

Accurate  tracking  of  the  KMAH  index,  k,  could  only  be 
accomplished  by  simultaneously  integrating  the  kinematic 
ray  tracing  equations  for  all  three  rays,  in  which  the  number 
of  sign  changes  of  the  distances  q,  are  accumulated  during 
the  time  step  of  the  integration.  In  this  approach  the 
computational  effort  is  roughly  the  same  as  the  procedure  of 
integrating  the  kinematic  system  once  and  the  vicinity  ray 
tracing  system  twice,  but  the  determination  of  quantities  qh 
q„  and  k  is  indirect  and  more  awkward.  Gajewski  & 
PSenCik  (1987)  have  applied  such  an  approach  to  determine 
dynamic  ray  tracing  quantities  in  anisotropic  media. 

In  this  application,  implementation  of  cither  the  standard 
dynamic  ray  tracing  equations  or  the  vicinity  ray  tracing 
equations  would  entail  the  evaluation  of  significantly  more 
complex  algebraic  expressions  at  each  time  step  in  the 
integration. 

3. 2. 5  Initial  conditions 

The  vicinity  rays  are  calculated  by  equations  (29),  while  the 
central  ray  is  calculated  by  the  ray  tracing  equations  in  the 
reference  coordinates  (e.g.,  Cartesian,  spherical,  and  etc.). 
Like  the  ray  tracing  equations,  equations  (29)  are  initial 
value  problems  in  non-linear,  first-order  differential 
equations,  which  can  be  solved  by  numerical  methods  such 
as  the  Runge-Kutta  method  (e.g.,  Forsythe,  Malcolm  & 
Moler  1977). 

Since  most  general  sources  can  be  constructed  from  a 
superposition  of  point-sources,  it  is  most  practical  to  specify 
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the  initial  conditions  of  the  VRT  system  for  a  point-source. 
These  are 

<7,1, -»  =  0  and  ij,U„=  if,1.  (31) 

where  superscript  /  denotes  the  initial  value  of  the 
parameter. 

In  a  superposition  of  Gaussian  beams,  there  are  two 
useful  choices  for  the  spacing  of  vicinity  rays.  Either  the 
vicinity  rays  can  be  chosen  to  be  mid-way  between  central 
rays  or  they  can  be  chosen  to  be  coincident  with  central 
rays.  In  the  mid-way  choice,  the  vicinity  rays  provide  a 
denser  sample  of  the  spatial  variation  in  traveitime  and 
hence  provide  a  more  accurate  estimate  of  the  phase  of  each 
complex  beam.  In  the  coincident  choice,  the  position  of  the 
estimated  vicinity  ray  can  be  compared  to  position  of  the 
exact,  kinematically  traced  ray.  thereby  checking  the 
accuracy  of  the  vicinity  ray  tracing  system  as  shown  in 
Section  3.4.  In  both  choices,  the  q[  are  independent  of  the 
Fresnel  beamwidth  described  in  Section  3.5. 


3.2.6  Relationship  between  VRT  and  DRT 

The  VRT  system  can  be  converted  into  the  DRT  system  by 
following  the  same  procedures  used  in  the  derivation  of  the 
DRT  system  (Cerveny  1985,  1987).  The  nght-hand  sides  of 
equations  (22)  can  be  divided  into  linear  terms  and 
non-linear  terms.  The  non-linear  terms  in  equations  (22) 
correspond  to  the  omitted  terms  in  DRT.  The  VRT  system 
in  equations  (22)  is  rewritten  as  follows  by  neglecting  higher 
than  the  second-order  terms  in  the  i-e,  plane  for  ».  /  ■  1,2, 
and  by  assuming  E  ■»  HV/v. 


dip, 

ds 


=  -Z**. 


(32) 


Equations  (32)  are  similar  to  the  paraxial  ray  tracing 
equations  (Cerveny  &  PSenfik  1979;  Cerveny  1987)  except 
for  the  scale  factor  A,.  VRT  includes  the  scale  factor  h, 
because  the  vicinity  (paraxial)  ray  is  evaluated  at  (j,  q,)  in 
VRT,  while  this  ray  is  evaluated  at  (s.  0. 0)  in  paraxial  ray 
tracing  and  in  DRT.  Equations  (32)  can  be  converted  into 
the  DRT  system  by  transforming  from  ray  centred  to  ray 
coordinates  (Cerveny  1985).  The  quantities  Q,t  and  P in  the 
DRT  are  obtained  from  q,  and  p,  in  the  VRT  for  i  and 
I  =  1,  2  by; 

( <?„  \  _  i  cos  i),  -sin  if,  w  q,\  _  /q,  cos  ij,\ 

'  Q,J  sin  i),  cos  t)(/'0 /  'q,  sini?,/ 

/cosij,  -sin  /p.cosija  (33) 

P,J  sin  r i,  cos  q,J '  0  /  p,  sin  q,  /  ’ 


where  i  *j.  i),  is  the  angular  difference  between  the  central 
and  the  vicinity  ray  in  the  i-e,  plane.  The  number  of 
equations  in  equations  (33)  is  doubled  compared  to 
equations  (32)  because  of  off-diagonal  terms.  The  physical 
meaning  of  q,  is  the  distance  from  the  central  ray  to  the 
vicinity  ray  measured  along  the  ray  centred  coordinate 
direction  <y  The  variation  of  q,  describes  the  change  in 
amplitude,  and  the  variation  of  q,  describes  the  geometry  of 
the  wavefront.  These  properties  of  q,  and  r\,  can  be  applied 
to  problems  such  as  two-point  ray  tracing,  estimation  of 
traveltime  near  a  central  ray,  and  calculation  of  geometric 
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spreading.  The  VRT  system  (29)  is  obtained  without  using 
any  paraxial  approximations.  Although  DRT  estimates  all 
rays  paraxially  close  to  a  reference  ray,  its  accuracy  rapidly 
deteriorates  as  the  distances  q,  increase  in  inhomogeneous 
medium.  VRT,  on  the  other  hand,  estimates  only  a  single 
ray  in  the  vicinity  of  a  reference  ray.  The  accuracy  of  this 
estimated  ray  remains  high  even  when  the  distant*  qt 
approach  the  scale  length  of  the  medium  (see  Section  3.4). 
In  either  traveltime  estimation  to  a  receiver  near  a  reference 
ray  or  seismogram  synthesis  by  superposition  of  Gaussian 
beams,  it  is  usually  possible  to  choose  a  spacing  of  vicinity 
rays  giving  a  higher  accuracy  in  estimated  traveltimes  than  is 
possible  using  paraxially  estimated  rays.  When  DRT  is 
sufficiently  accurate  in  a  slowly  varying  medium ,  equations 
(33)  can  be  used  to  convert  between  VRT  and  DRT. 

Because  the  velocity  of  the  vicinity  rays  are  not  expanded 
at  the  point  of  the  central  ray  (,v,  0. 0)  in  the  VRT,  the  VRT 
system  requires  three  velocities  i>,  V,  and  V,  and  their 
derivatives  at  points  (s,  0,  0),  (s.q,.  0)  and  (s,  0,  q2).  This 
requirement  may  increase  computing  time  compared  to  the 
DRT  system,  which  requ  is  only  one  velocity  v  and  its 
second  spatial  derivatives  at  point  (j,  0, 0).  The  VRT 
system,  however,  does  not  require  the  computation  of 
second  spatial  derivatives  and  avoids  approximations  in 
calculating  the  velocities  of  the  medium  along  the  vicinity 
rays.  Note  that  the  equations  for  q,  and  q2  depend  on  the 
velocity  of  the  medium  along  the  vicinity  ray,  V,  or  V2  rather 
than  on  the  velocity  of  the  medium  along  the  central  ray, 
V(r,  0, 0)  =  v.  For  a  velocity  model  specified  in  fixed 
Cartesian  (or  spherical)  coordinates,  the  velocity  V,  and  its 
derivatives  V,V(  can  be  calculated  by  transforming  the 
positions  of  the  vicinity  ray  in  ray  centred  coordi>:~.:s 
(0,  qit  0)  and  (0,  0,  q2)  to  fixed  Cartesian  (or  spherical) 
coordinates. 

Because  the  VRT  system  calculates  q,  and  q,  values  by 
using  V,  and  V,%  it  is  not  necessary  to  employ  the  method  of 
matching  paraxial  phase  (Cerveny  &  Hron  1980)  to 
determine  new  initial  conditions  on  q,  and  q,  when  vicinity 
rays  are  transmitted  through  or  refiected  by  discontinuities. 
Since  second  spatial  derivatives  of  velocity  are  not  used  in 
the  VRT  system,  no  jumps  in  q ,  or  q,  are  induced  by 
discontinuities  in  velocity  gradient.  At  first-order  discon¬ 
tinuities  in  velocity,  new  initial  conditions  on  q,  and  q,  are 
computed  by  simply  applying  Snell's  law  to  both  the  central 
ray  and  the  vicinity  ray. 

The  differences  between  VRT  and  DRT  can  be 
summarized  as  follows. 

(1)  Wavefront  curvature,  geometric  spreading,  traveltime 
to  a  point  near  a  central  ray,  and  an  iterative  solution  to  the 
two-point  ray  tracing  problem  can  be  obtained  using  either 
the  DRT  or  VRT  systems. 

(2)  VRT.  if  desired,  can  be  converted  into  DRT  in 
regions  of  a  medium  where 


Likewise,  DRT  can  be  converted  into  VRT  in  more  rapidly 
varying  regions  of  a  medium  using  equations  (33)  and 
choosing  an  initial  q,  to  continue  the  integration  using  VRT, 

(3)  Since  both  VRT  and  DRT  require  ray  centred 
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coordinates,  S- wave  polarization  is  handled  in  VRT  using 
the  same  procedure  as  in  DRT  (terveny  1987). 

(4)  The  geometric  spreading  of  a  wavefront  in  VRT  is 
determined  from  the  conservation  of  energv  flu*.  |n  DRT  it 
is  determined  from  an  approximate  solution  of  the  transport 
equation. 

(5)  DRT  requires  specification  of  first  and  second  spatial 
derivatives  of  the  medium  and  jump  conditions  at 
discontinuities  in  velocity  gradient.  VRT  requires  specifica¬ 
tion  of  only  the  first  spatial  gradient  of  the  medium  and  does 
not  require  jump  conditions  at  gradient  discontinuities. 

(6)  DRT  can  be  used  to  estimate  the  trajectory  and 
traveltime  of  all  rays  paraxially  close  to  a  reference  ray. 
VRT  can  be  used  to  estimate  the  trajectory  and  traveltime 
of  a  single  ray  in  the  ray  centred  coordinate  system  of  a 
reference  ray.  This  estimate  is  much  more  accurate  than  that 
possible  with  DRT  and  remains  accurate  even  as  the 
distances  q,  approach  the  scale  length  of  the  medium. 


3.3  Computation  of  traveHime  near  a  central  ray 

The  computation  of  the  traveltime  to  a  receiver  near  a 
central  ray  is  just  as  simple  in  the  VRT  system  as  in  the 
standard  method  of  DRT  using  the  paraxial  approximation. 
Fig.  3  illustrates  the  calculation  of  the  traveltime, 
t(s,  nt,n2)  at  point  N(s,nun2).  The  determination  of  s 
and  n ,  for  a  specified  point  N  in  ray  centred  coordinates  is 
important  in  obtaining  accurate  estimates  of  traveltime  and 
amplitude  of  the  vicinity  ray  with  respect  to  a  central  ray. 
The  rough  approximations  contained  in  the  standard 
paraxial  technique  may  produce  spurious  oscillations  in  the 
superposition  of  Gaussian  beams  (e.g.,  Madariaga  1984)  and 
break  down  if  the  central  ray  is  far  from  the  receiver.  Here, 
we  describe  an  alternative  method  for  the  determination  of 
a  specified  point  in  ray  centred  coordinates.  We  begin  by 
writing  the  traveltime  field,  r(r, /i|,n2)  of  the  specified 
point  N  (e.g.,  receiver),  in  the  ray  centred  coordinates  as 

r(r,  n , ,  /i2) =  r(s )  +  A  r.  (34) 

Ar  is  the  traveltime  difference  direction  between  the  points 
S  and  N  (Fig.  3).  The  traveltime  difference  Ar(x,  nt,n2) 
between  S  and  N,  is  obtained  from 


Ar  = 


(35) 


where  An  is  the  distance  from  N  to  a  point  D  on  the 
wavefront  measured  along  a  normal  to  the  wavefront,  and  V 
is  the  velocity  of  the  vicinity  ray  at  point  (s,  n2).  e  =  ±1 

or  0  depending  on  the  shape  of  the  wavefront  along  the 
direction  whose  projected  line  passes  the  points  (s,  0,  0)  and 
( s ,  n,,  n2)  in  the  e,-t\  plane.  Let  us  assume  that  the  radius 
R,  (or  curvature  K ,)  of  the  wavefront  along  the  e,-axis  does 
not  change  in  the  neighbourhood  of  the  point  S(s,  0,  0). 
When  the  point  is  located  along  the  axis,  e,  or  e2,  the 
distance  An,  is  simply  approximated  as  shown  in  Fig.  3,  and 
is  given  by 

An,  -  VRjT^  -  R,.  (36) 

The  radius  R,  at  the  point  (s,  n,.  0)  or  (r,  0,  n2)  is  given  by 
equation  (24).  The  error  in  equation  (36)  depends  on  the 
variation  of  R,  along  the  e,-axis.  Let  q\  denote  the  normal 


distance  to  the  central  ray  from  point  B,  where  B  is  the 
intersection  of  the  wavefront  of  the  central  ray  and  vicinity 
ray  (Fig.  3).  The  curvature  (K,)  and  the  radius  of  curvature 
of  the  wavefront  (R,)  of  the  principal  axis,  i„  at  point 
(r,  0, 0)  can  be  expressed  as  follows  (Cerven^  1981)  by  using 
equation  (30): 

A  _  _L  —  jL  _  Si  —  /j7) 

‘  ~K,  V(p,  ~  sin  q,  tanij/ 

Let  ifi  denote  the  distance  from  S  to  B  along  the  wavefront. 
q,  is  the  distance  along  the  wavefront  corresponding  to  q,  in 
the  ray  centred  coordinates.  The  wavefront  coordinates 
(r,  qu  q2)  are  defined  in  Fig.  3.  The  relation  between  q,  and 
q,  can  be  represented  using  R,: 


tan  q,  q, 


(38) 


Equation  (38)  determines  the  Jacobian  1  between  the  ray 
centred  coordinates  and  the  wavefront  coordinates: 


<?I<?2=-^1<?2>  (39) 

where 


J  = 


9l»h 

tan  tan  q2 


Equation  (38)  shows  that  the  curvature  of  a  wavefront  K,  or 
the  radius  of  curvature  of  the  wavefront  R,  can  be  written  as 
a  simple  function  of  q ,  and  q,.  When  the  vicinity  is  located 
on  the  e,-axis,  At  is  simply  obtained  by  substituting 
equation  (36)  into  (35): 


Ar=  e, 


An, 

V, 


VRf  +  nf-R, 
K 


(40) 


where 


e,  =  1  for  q,  x  q,  >0:  convex  wavefront, 

e,  =  0  for  q,xq,=  0:  planar  wavefront, 

f,  =  - 1  for  q,  x  rj,  <  0:  concave  wavefront. 

To  facilitate  comparison  with  standard  DRT  in  2-D, 
equation  (40)  can  be  expanded  as  follows  for  n/R  «  1: 


The  traveltime  of  the  vicinity  ray  r(s,  n)  is  approximated  as 
follows  by  substituting  equation  (41)  into  (34): 


r(.r,  n)«  r(j)+  \( 


VR 


(42) 


Equation  (42)  is  similar  to  that  in  DRT.  The  factor  l/V'ff  in 
equation  (41)  corresponds  to  M  (equation  9)  in  the 
expression  for  standard  DRT: 


1  sin  q  _  p  p 

VR  Vq'  q'  q-Ans\nq 


M  +  M  An  sin  q. 


(43) 


The  last  term  of  the  right-hand  side  in  equation  (43)  is 
neglected  in  the  paraxial  approximations  of  DRT.  q'  is  the 
normal  distance  between  C  and  D  in  Fig.  3.  As  shown  in 
Table  3  in  Section  3.4,  equation  (40)  gives  more  accurate 
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computational  results  than  the  paraxial  approximation  in 
DRT  does.  When  the  vicinity  ray  is  located  at  the  general 
point.  S'(s.  nt,n2)  in  the  ex-e2  plane.  Ar  can  be  obtained 
by  calculating  the  radius  of  curvature  of  the  wavefront.  R,  at 
the  point  N(s.  nx,  n2)  along  the  new  principal  axis  e.  The 
angle  p  is  an  angular  difference  between  the  e,  and  e  axes  in 
the  i\-f2  plane  (Fig.  4).  By  using  the  new  principal  axis  e. 
the  point  N(s,  n,.  n2)  in  3-D  can  be  described  as  N(s,  n)  in 
2-D.  The  distance,  n,  between  the  points  S  and  N  in  the 
e,-e,  plane  is  given  as 

n  =/»,  cosp  +  /t,sin/r,  (44) 

where 


The  radius  of  curvature  of  the  wavefront,  R.  at  point  N(s,  0) 
in  the  i-e  plane  is  obtained  by  using  E  vr’s  theorem 
(Cerveny  &  Ravindra  1971): 


11,  1  .  , 

-  =  — cos-p+  — s.n-p. 


(45) 


Finally,  Ar  at  the  point  N(s, «,)  =  N(s,  n)  can  be 
calculated  by  using  equations  (43),  (44)  and  (45)  in  the  i-e 
plane. 


VF+Tr-R  R 

4r-‘ - V - •  wh're  <  =  Ri- 


(«) 


r(r.  n,)  in  equation  (34)  can  be  described  in  terms  of  known 
point  j  =  j„  along  the  central  ray.  The  quantity  r(s)  can  be 


Figure  4.  New  principal  axis  e  and  the  i-e  plane:  The  new  principal 
e  axis  is  constructed  by  rotating  the  e,-axis  through  the  angle  p 
which  is  described  in  terms  of  n,  and  n>  The  point  N(s,  n,,  n2)  in 
f-D  can  be  described  as  N(s,n)  in  2-D  by  constructing  the  i-e 
plane.  The  distance  n  in  the  i-e  plane  is  obtained  by  using  equation 
(51). 


expanded  with  respect  to  r(j„)  by  using  a  Taylor  expansion 
about  r„.  Terms  higher  than  second  order  are  negligible  and 
will  be  neglected 


r(r)  =  r(s„)  + 


3r(r) 

3s 


«~IO 


1  d2X  (s) 

+  2  <?r 


(r-j„)2  +  ---. 


(47) 


It  is  easy  to  see  that 


3r(i) 


ds 


1 

v(s„) 


and 


a2x(s) 


ds2 


vj*0 ) 
v2(su) ' 


(48) 


Combining  the  expressions  (34),  (46),  (47)  and  (48),  the 
traveltime  t(r, «,)  is  approximated  by 


x(s,  tij)  =  t(j0)  +  (s  ~s„)  -  (s  ~  •'o)2  +  Ar.  (49) 

x(s,  n,)  indicates  the  traveltime  of  a  specified  point  N,  such 
as  a  receiver  point,  in  ray  centred  coordinates.  Note  that 
although  a  Taylor  expansion  has  been  used,  it  is  a  Taylor 
expansion  along  the  direction  of  the  central  ray  rather  than 
along  a  direction  perpendicular  to  the  central  ray.  The 
standard  Gaussian  beam  and  paraxial  ray  methods  make  a 
Taylor  expansion  of  traveltime  in  the  direction  perpendicu¬ 
lar  to  the  central  ray  as  well.  In  equation  (49)  it  is  usually 
possible  to  select  close  to  s,  minimizing  the  errors  in 
making  a  Taylor  expansion  along  the  central  ray. 

Two  important  modifications  from  standard  DRT  will 
make  the  estimated  traveltime  more  accurate  in  VRT.  First, 
the  radical  in  equation  (46)  is  not  expanded.  Second,  the 
wavefront  curvature  is  calculated  at  the  vicinity  ray  rather 
than  the  central  ray.  By  judicious  choice  of  the  spacing  of 
central  rays  and  associated  vicinity  rays,  it  is  possible  to 
more  accurately  describe  the  true  shape  of  the  wavefront. 
The  improved  accuracy  possible  with  the  VRT  system  has 
been  demonstrated  in  several  tests. 

Estimation  of  the  traveltime  for  a  vicinity  ray  becomes 
less  accurate  as  the  distance  is  increased.  The  error,  for  a 
specified  point  (e.g.,  receiver),  can  be  reduced  by  estimating 
the  traveltime  from  the  closest  ray  to  the  receiver,  normal  to 
the  wavefront  passing  through  the  receiver.  In  the  case  of 
multipathing,  several  wavefronts  may  pass  through  a 
receiver.  Each  wavefront  is  classified  by  its  KMAH  index 
(see  Section  3.5)  and  the  signs  of  qt  and  r),  of  the  rays 
belonging  to  the  wavefront.  A  traveltime  is  then  estimated 
for  each  wavefront  passing  through  the  receiver. 


3.4  The  accuracy  of  the  VRT  and  DRT  systems 

The  accuracy  of  the  DRT  and  VRT  systems  was  tested  in 
computations  using  the  structure  shown  in  Fig.  5.  In  this 
model,  the  velocity  gradient  is  discontinuous  at  zlt  =  - 10  km 
but  velocity  itself  remains  continuous.  Because  of  a 
discontinuity  in  the  curvature  of  wavefront  at  zd,  the  scale 
factor  h  in  ray  centred  coordinates  is  also  discontinuous  at 
z,,.  The  scale  factor  h  at  a  point  (*,,  z)  is  calculated  as 
follows  when  the  central  ray  is  below/above  zd,  and  its 
vicinity  ray  is  above/below  zlt : 

h  =  l  +  Wl^*+W2^, 
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Flftit  5.  A  laterally  homogeneous  acoustic  model  with  density 
constant  in  depth:  the  gradient  of  the  velocity  has  a  discontinuity  at 
2  =  10  km.  The  velocity  of  the  medium  is  given  by  v(z)  = 
2.5  +  0.1  x|z|  for  za-lOkm,  and  u(z)  =  3.5  + 0.4  x(|z|-  10)  for 
z  s  -10km. 

where 


*2  -  2,1 

w  — 

Zi  -  z,r 

Zj-Z, 

t  ns  — 

*2-Z| 

dq 

z2  =  zt-q  sin  or.  sin  a  =  — . 


v  is  the  velocity  of  the  central  ray  at  the  point  S,  and  V  is 
the  velocity  of  its  vicinity  ray. 

The  quantities  of  vicinity  rays,  q  and  p,  and  estimates  of 
traveltime  were  calculated  using  equations  (29)  and  (40)  as  a 
test  of  the  VRT  system,  and  using  equation  (30)  and  (2)  as  a 
test  of  the  DRT  system.  In  DRT,  the  phase  matching 
method  (terven^  1981)  is  applied  at  the  discontinuity  in 
velocity  gradient,  zd  =  -10  km,  while  it  is  not  required  in 
the  VRT.  Figs  6  and  7  show  the  central  ray  Qq,  and  its  true 
vicinity  rays,  Q,  for  the  initial  angular  difference  » j°  =  «, 
where  i=  1,2,3.  r f=i  means  that  the  vicinity  ray’s  take-off 
angle  is  i  degrees  away  with  respect  to  a  particular  central 
ray.  Therefore,  if  there  are  no  errors  in  VRT  and  in  DRT, 
the  vicinity  ray  path  must  be  that  of  the  prior  (or 
subsequent)  central  ray.  In  Figs  6  and  7,  solid  lines  indicate 
the  central  ray  path  and  its  true  vicinity  ray  paths,  while 
dotted  lines  represent  its  estimated  vicinity  ray  paths  from 
VRT,  Q.'f,  and  from  DRT,  Qf.  The  true  vicinity  ray  paths 
are  calculated  from  kinematic  ray  tracing,  by  using  a 
take-off  angle  or,,  + 1,  where  a-,,  is  the  take-off  angle  of  the 
central  ray  at  the  source.  Therefore,  the  difference  between 
f2„  solid  line,  and  Q.'f  (or  Qf),  dotted  line,  represents 
errors  in  VRT  (or  in  DRT).  As  shown  in  Figs  6  and  7,  the 
errors  of  q  and  p  in  VRT  are  considerably  smaller  than 
those  in  DRT.  Tables  1  and  2  give  more  detailed 
information,  comparing  true  values  of  q  and  p  with  the 
estimated  values  from  the  VRT  and  DRT.  The  symbol  A 
denotes  the  difference  between  the  true  value  and  its 
estimated  value. 

As  shown  in  Tables  1  and  2,  the  accuracy  of  VRT  is 
greater  than  that  of  DRT.  This  is  because  VRT  does  not 
make  any  paraxial  approximations  and,  hence  should 
remain  accurate  for  much  greate  distances  q,  from  the 
central  ray.  The  errors  caused  by  approximation  of  the  scale 
factor  h  in  VRT.  however,  are  inescapable.  Errors  will  grow 
in  VRT  as  q,  increases  but  at  a  much  slower  rate  than  in 
DRT.  The  trajectory  and  traveltime  of  vicinity  rays  can  be 
determined  exactly  by  kinematic  ray  tracing  with  geometric 


True  vicinity  ray  paths  and  their  estimated  ray  paths  from  the  VRT:  the  true  vicinity  ray  paths  (solid  lines)  for  =  /'  =  1,  2  and  3 
degrees  (Q,)  and  the  estimated  vicinity  ray  paths  (dotted  lines)  (Q,v)  are  plotted  to  test  the  accuracy  of  q  and  p  in  the  VRT  system.  The 
differences  between  the  solid  lines  and  the  dotted  lines  represent  the  errors  in  the  VRT 
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Figwc  7.  True  vicinity  ray  paths  and  their  estimated  ray  paths  from  the  DRT:  the  true  vicinity  ray  paths  (solid  lines)  for  tj",  f  =  I,  2  and  3 
degrees  (Q,)  and  the  estimated  vicinity  ray  paths  (dotted  lines)  (Q,°)  are  plotted  to  test  the  accuracy  of  q  and  p  in  the  DRT  system.  The 
differences  between  the  solid  lines  and  the  dotted  lines  represent  the  errors  in  the  DRT. 


Table  1.  Estimated  values,  qv  and  pv .  from  the  VRT. 


1° 

True  p 

True  q 

PV 

< 

1V 

Aqv 

I* 

-3  995E-3 

16  6820 

-3.392E-3 

-0  653E-3 

16.5673 

0.1147 

2* 

-8.58IE-3 

29.0914 

-6.635E-3 

•  I.946E-3 

28.3819 

0.6095 

3* 

-1418E-2 

38.6256 

-1.032E-2 

-3  860E-3 

36.7859 

0.8397 

Table  2.  Estimated  values,  q"  and  p“,  from  the  DRT. 
T)‘  Tni  ep  Trueq  p"  A  p"  q" 

&q" 

1*  -3  995E-3 

16.6820 

-2.362E-2  1  963E-2 

18  2923 

-1.6103 

2*  -8.571E-3 

29.0914 

-4  723E-2  3  865C-2 

36.5790 

-7  4876 

3*  -I418E-2 

38.6256 

-7  082E-2  5  564E-2 

54.8545 

-16.2289 

spreading,  wavefront  curvature,  and  KMAH  indices 
calculated  as  described  by  Gajewski  &  PSenCik  (1987).  In 
this  procedure,  the  trajectory  of  vicinity  rays  is  error  free, 
but  since  the  spreading  and  curvature  are  calculated  from 
differencing  closely  spaced  rays,  they  potentially  have 
greater  error  than  those  obtained  from  the  numerically  more 
stable  VRT  procedure. 

Table  3  shows  the  traveltime  estimation  of  the  vicinity  ray 
from  VRT  and  from  DRT,  with  respect  to  the  traveltime  of 
the  central  ray.  As  shown  in  Fig.  8,  the  true  travcltimes  of 
the  vicinity  rays  are  obtained  from  kinematic  ray  tracing.  To 
avoid  errors  generated  along  the  central  ray  path,  the 


Flfart  8.  Traveltime  estimation  from  the  VRT  and  from  the  DRT:  the  true  traveltimes  of  the  vicinity  rays  are  obtained  from  the  kinematic  ray 
tracing.  The  traveltimes  of  the  vicinity  rays  at  /? , ,  R2  and  Rx  are  estimated  at  the  point  St,  S2  and  S2  on  the  central  ray  to  avoid  errors 
generated  along  the  central  ray  path.  The  true  q  and  p  are  used  in  the  traveltime  estimation  for  both  the  VRT  and  the  DRT. 
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Table  3.  Traveltimc  estimation  from  the  VRT  and  the  DRT 


n' 

t  rue  T 

Estimated  Tl 

AT1 

Estimated  f" 

A7" 

r 

19 1(725 

19  0701 

0.0024 

19  (.382 

0  0343 

2* 

1X41% 

18  4027 

0  0169 

18  2832 

0  1364 

V 

17X480 

17  7886 

0  0594 

17  5267 

0  3213 

traveltimes  of  vicinity  rays  at  Rlt  R2  and  are  estimated  at 
points  S,.  S2  and  on  the  central  ray  path.  The  true  values 
of  q  and  p  are  used  in  the  traveltime  estimation  of  vicinity 
rays.  The  traveltime  estimation  of  the  vicinity  ray  from  VRT 
is  calculated  by  using  equation  (29),  and  that  from  DRT  is 
calculated  by  using  equation  (2).  As  shown  in  Table  3,  the 
estimated  traveltimes  from  VRT  are  more  accurate  than 
those  from  DRT.  When  the  estimated  q  and  p  are  used  in 
the  traveltime  estimation,  the  estimated  traveltime  errors  in 
DRT  greatly  increase  compared  to  those  estimated  from 
VRT  due  to  larger  errors  in  q  and  p  in  DRT. 

3.5  The  synthesis  of  seismograms  by  superposition  of 
Gaussian  beams 

The  zeroth-order  high-frequency  asymptotic  solution  to  the 
wave  equation  in  generalized  coordinates  is  given  in 
equation  (12).  Let  us  consider  a  point  S,  located  close  to  the 
central  ray,  specified  by  the  ray  centred  coordinates, 
(s.  n,,n2),  that  is,  S  =  (s,nx,n2).  The  zeroth-order 
asymptotic  solution  to  the  reduced  wave  equation  in  the  ray 
centred  coordinates  can  be  expressed  in  the  form 

g(S,  to)  =  A(S)e~'""<s,(-i)k  sgn  (to).  (50) 

The  amplitude  function  A  is  complex  and  can  be  determined 
by  applying  the  conservation  of  energy  flux  and  a 
normalization  condition  (Beisei  1969,  p.  156;  Gasiorowicz 
1974,  p.  45).  r  is  the  traveltime  of  the  central  ray.  k  is  the 
value  of  the  KMAH  index  whose  value  is  increased  by  one 
whenever  the  sign  of  q,  changes  along  the  ray.  The  KMAH 
index  represents  the  x/2  phase  shift  whenever  the  ray 
touches  a  x-caustic  (</ ,  =  0  or  q2  =  0)  (Chapman  & 
Drummond  1982). 

3.5. 1  Source-time  functions 

For  a  source-time  function  f(t)  specified  as  the  real  part  of 
an  analytic  function  v(f ).  the  wavelield  is  given  by 
evaluating  a  convolution: 

u(S,  t)  =  &  |g(.S\  /)  *  y(0l 

=  g(S,  <o)y(to)e  ,uu  dw^j  n  (51) 

and  substituting  equation  (50)  for  g(S,  to)  gives 

u(S,  T)  -  A-  f  (-/)*  sgn  (u))y(w)e  "dto 
n  J  , 

Jfc.|y(r-r)(-i)*|.  (52) 

71 

v(i)  is  the  analytic  time  scries  represented  bv 

y(t)  -/(f)  -  th{t),  (S3) 

where  /(/)  and  /i(f)  are  a  Hilbert  transform  pair  (Bracewell 
1978)  The  analytic  function  corresponding  to  anv  realistic 


source-time  function  can  be  constructed  by  choosing  y(0  to 
be  a  generalized  delta  function  and  convolving  that  function 
with  a  particular  source  time  function.  Some  possible  forms 
for  y(i)  are  for  /(/)  equal  to: 


(1)  a  delta  function  (Chapman  1978;  Chapman  & 
Drummond  1982): 

y(t)  =  &(t)-i  — :  (54) 

Jit 

(2)  a  Gaussian  wavelet  (Cerveny  1983): 


y(0~ 


<u  x>Ytl*£-t'oiiYtili  -  no! 

Wry 


(55) 


where  y  controls  the  width  of  a  Gaussian  envelope  with 
respect  to  the  prevailing  angular  frequency  w:  and 

(3)  a  resonance  function  (Madariaga  &  Papadimitriou 
1985): 


At 

tJ  +  (At)2 


r  +  (At)2 


(56) 


where  At  is  the  sampling  interval  for  the  discrete  time  series 
>•(')• 

The  Gaussian  wavelet  equation  (55)  is  useful  for  simulating 
a  narrow-band  source,  while  equation  (56)  is  useful  for 
simulating  broad-band  responses.  Equations  (54),  (55),  (56) 
can  be  constructed  to  be  a  generalized  delta  function  by 
requiring 

f»y(t)<ft-l.  (57) 


3.5.2  Beamwtdlhs  and  the  beam  Fresnel  volume 

In  beam  methods,  every  ray  can  contribute  to  the  wavefield 
at  the  receiver,  but  the  contributions  are  dependent  on  the 
particular  beamwidths.  An  important  problem  in  the 
superposition  of  beams  is  the  determination  of  the 
beamwidths.  The  Fresnel  volume  is  one  possible  form  of  a 
physical  bcamwidth,  because  the  principal  contribution  to 
the  amplitude  and  the  phase  at  the  receiver  is  contained  in 
the  tays  within  the  beam  Fresnel  volume.  Fresnel  diffraction 
and  the  Fresnel  volume  are  defined  in  optics  (Jenkins  & 
White  1937;  Stone  1963)  and  in  geophysics  (Kravtsov  & 
Orlov  1980;  Cervenv  1987)  to  explain  diffraction.  An 
approximate  Fresnel  volume  can  be  defined  in  the  beam 
method,  which  we  refer  to  as  the  ‘beam  Fresnel  volume’.  A 
vicinity  ray  contained  on  the  surface  of  the  beam  Fresnel 
volume  is  half-period  ahead  (or  behind)  in  traveltime  to  that 
of  the  central  ray.  For  a  given  frequency,  the  cross-section 
of  the  beam  Fresnel  volume  perpendicular  to  the  central  ray 
at  a  point  s  can  be  estimated  from  all  vicinity  rays  at  s  a 
half-period  ahead  or  behind  the  traveltime  of  the  central 
ray.  The  beam  Fresnel  volume  is  not  the  same  as  the  Fresnel 
volume  defined  by  Kravtsov  &  Orlov  (1980),  but  its 
cross-section  closely  approximates  the  cross-section  of  the 
exact  Fresnel  volume  along  most  of  the  rav  path  of  the 
minimum  time  ray. 

The  tiaveltime  difference  between  the  central  ray  and  a 
vicinity  rav  on  the  surface  of  the  beam  Fresnel  volume  in  the 
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i-e,  plane  is  given  as 
A  r  =  r(v.  T, ,  0)  -  r (s,  0,  0)  =  y, 

(58) 

Ar  =  r(s,  0,  F2)  -  t(j,  0,  0)  =  y, 

where  y  is  the  half-period.  Equation  (58)  states  that  a  point 
on  the  surface  of  the  Fresnel  volume  is,  F, .  0)  or  is,  0,  F:) 
has  a  half-period  time  difference  with  respect  to  the  point 
(s.  0.  0)  on  the  central  ray.  Formulae  for  the  beamwidths  F, 
along  the  e,  axes  arc  obtained  by  substituting  equation  (40) 
into  (58). 

F\  =  Jy*V\  +  2  yVx  =  VyJF;  +  2yK,/?t. 

V  tan  t}  i 

I -  (59) 

F,  =•  V  fV\  +  2  yV2  — -  =  Vy^F;  +  2yV,R2. 

V  tan  i)  2 

At  a  fixed  point  q,  and  rj,  may  vary  according  to  the  initial 
conditions  chosen  for  vicinity  ray  tracing.  The  ratio 
q,/tan  r)„  however,  is  nearly  independent  of  the  initial 
conditions.  The  definitions  (59)  ensure  that  the  beam  widths 
are  always  equal  to  the  width  of  the  approximate 
cross-section  of  the  Fresnel  zone  for  a  particular  frequency, 
independent  of  the  initial  conditions  chosen  for  vicinity  ray 
tracing.  For  high  frequency  (small  y).  the  beamwidth  given 
by  (59)  is  approximately  proportional  to  \J2 yR,V,  = 

Equation  (59)  is  the  same  as  the  classical  definition  of 
Fresnel’s  half-period  zones  (e.g.,  Jenkins  &  White  1937). 

Using  Fresnel  beamwidths,  the  amplitude  function  A  at  a 
specified  point  N  will  be  described  as  a  generalized  Gaussian 
function  of  the  form 


A(5’)  =  C'exp 


where  the  beamwidth  F,  is  the  half-width  of  the  beam 
Fresnel  volume.  With  this  amplitude  distribution,  energy 
along  the  beam  Fresnel  volume  of  half-width  F,  is 
proportional  to  1/e'  and  its  amplitude  is  proportional  to  1/e 
with  respect  to  the  central  ray.  If  the  expression  for  the 
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Figart  9,  The  vicinity  ray  and  the  beam  Fresnel  volume:  the 
distance  q,  between  the  central  and  the  vicinity  ray  is  generally 
smaller  than  the  beamwidth  F,  from  the  central  ray  to  the  beam 
Fresnel  volume.  The  parameter  at  point  A,  and  q,  =0  at 
point  H.  which  correspond  to  the  v  and  x  caustics  (Chapman  & 
Drummond  1982). 


amplitude  function  is  viewed  as  a  generalized  delta  function. 
C  can  be  chosen  from  a  normalization  condition  in  space 
and  time  including  the  reflection  and  the  transmission 
coefficients.  This  will  guarantee  that  A[S)  will  satisfy  the 
transport  equation  and  that  equation  (52)  will  satisfy  the 
high-frequency  equation  of  motion, 

Generally.  F,  is  chosen  to  be  equal  to  the  half-width  of  the 
beam  Fresnel  volume  as  defined  in  equations  (59).  A 
modification,  however,  is  necessary  near  y  caustics  (Fig.  9), 
where  x  and  y  caustics  are  as  defined  in  Chapman  & 
Drummond  (1982).  Note  that  the  formulae  (59)  are  regular 
at  x  caustics,  where  q,  and/or  q2  vanish.  At  y  caustics, 
where  q,  and/or  rj2  vanish,  the  formulae  (59)  go  to  infinity. 
The  y  caustics  correspond  to  plane  waves.  Since  some 
Gaussian  windowing  is  always  desirable  to  exclude 
contributions  from  vicinity  rays  having  errors  at  large 
distances  q„  we  take  the  beamwidth  in  the  vicinity  of  a  y 
caustic  to  be  equal  to  the  mean  Fresnel  beamwidth  in 
regions  adjacent  to  the  y  caustic.  This  modification  near  y 
caustics  does  not  violate  energy  conservation  and  the 
normalization  condition. 


3.5.3  Energy  conservation  and  normalization  of  beams 

Using  equations  (52)  and  (60)  the  complex  displacement  u 
of  a  beam  specified  at  a  point  N  in  the  ray  centred 
coordinates  is 

u(S,  t)  -  Cexp  [-(^)2-  (^)2]  -  r)l.  (61) 

where  the  amplitude  factor  C  is  obtained  by  using  the  law  of 
conservation  of  energy  flux  and  a  normalization  condition. 
This  approach  for  determining  C  differs  from  the  approach 
followed  in  the  standard  Gaussian  beam  technique,  where  C' 
is  obtained  by  evaluating  the  superposition  integral  by 
stationary  phase  and  requiring  the  result  to  reproduce  ray 
theory  in  regions  where  it  is  valid.  In  contrast  to  the 
standard  Gaussian  beam  method,  beams  are  interpreted  as 
the  probability  of  finding  a  ray  at  a  given  point  and  time. 
This  probability  distribution  is  assumed  to  be  a  Gaussian 
distribution  whose  unit  area  is  always  1  with  respect  to  n, 
and  t.  This  constraint  guarantees  that  the  energy  of  a  beam 
(ray  tube)  is  conserved  with  respect  to  space  and  time. 
Conservation  of  energy  flux  in  the  VRT  system  is  identical 
to  the  requirement  that  the  determinant  of  the  propagator 
matrix  in  the  DRT  system  is  constant  along  the  ray 
(Liouville’s  theorem)  (terveny  &  PSenCik  1983;  Kim  1986; 
Klimeg  1988).  The  wave  function  u(S,  t)  then  describes  the 
probability  of  finding  a  ray  with  a  statistical  state,  which  is 
characterized  by  u.  Because  the  total  energy  in  a  beam  is 
conserved  along  the  wavefront,  it  is  necessary  to  transform 
from  the  ray  centred  coordinates  to  the  wavefront 
coordinates.  The  Jacobian  between  the  ray  centred  and  the 
wavefront  coordinates  is  given  in  equation  (39).  The 
conservation  law  of  energy  flux  and  the  normalization 
condition  imply  that  the  probability  (P)  of  finding  a  ray 
within  a  given  space  is 


P(S,  0  =  (u  I  “)i.o,o 

Ooe. 

dqxdq2Ju(S,  t)u*(S,  t)  =  1. 

-eo 
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where  (he  symbol  *  denotes  the  complex  conjugate  and  J  is 
the  Jacobian  between  the  ray  centred  and  the  wavefront 
coordinates  (equation  39)  The  probability  of  finding  a  ray  is 
a  similar  concept  to  the  density  of  rays.  If  the  integration 
limits  in  equation  (62)  are  taken  over  a  fixed  volume  of 
space,  then  a  focusing  tegion  would  give  a  large  P  value, 
while  a  defocusing  region  would  give  a  small  P  value.  The 
constant  C  in  equation  (62)  is  determined  by  solving 
equation  (62)  and  by  considering  the  normalization  factor  D 
for  a  radiation  pattern  at  the  source  and  <1>  for  a  product  of 
reflection  and  transmission  coefficients.  The  integral  of 
equation  (62)  yields 

P(i',/)  =  C'  =  (63) 


The  constant  D  depends  on  the  take-off  angle  6  and  azimuth 
P  (Aki  &  Richards  1980.  p.  82;  KlimeS  1984).  Reflection 
and  transmission  coefficients  are  introduced  by  many 
authors  (Aki  &  Richards  1980;  terveny  1985,  etc.).  The 
product  d>  is  generally  a  complex  value,  and  is  given  as 

=  n  / 


where  L,  i?  a  reflection/transmission  coefficient  at  an  tth 
interface.  The  constant  C  is  then  obtained  as 

C-'f20’*’ 


yj.th]  F2  represents  the  geometrical  spreading  of  the 
wavefront.  A(S)  can  be  rewritten  as  follows  by  substituting 
the  expres:  ion  for  C  into  equation  (60); 


,  ,  \f20t>  { 

^  s-TTpe!lP 
W/tfj  F.  t 


(64) 


The  calculation  of  the  traveltime  r(r.  n,)  in  (50)  is  given  in 
equations  (49).  The  final  wavcheld  of  a  ray  at  a  specified 
point  is  then  obtained  as 


u(.S.  l)  •= 


V2  D 

F, 

x  exp  [-(-•)  -  (j?)  ]  |«W-/)\v(f  -  r)|.  (65) 


The  form  of  equation  (05)  can  be  shown  to  be  quite  similar 
to  the  expression  for  a  standard  Gaussian  beam  when 
paraxial  approximations  are  substituted  for  the  phase  and 
the  expression  of  the  half-width  of  the  paraxial  volume  is 
substituted  for  h]  The  approach  and  concepts  used  in 
dcuvitig  the  vicinity  ray  tracing  system,  however,  are  quite 
different  from  those  used  in  the  standard  Gaussian  beam 
method.  The  number  of  equations  required  m  this  method  is 
nine,  by  contrast  to  21  in  the  standard  Gaussian  beam 
method.  This  method  uses  exact  positions  of  vicinity  rays 
while  the  standard  Gaussian  beam  method  uses  estimated 
values  based  on  a  Taylor  expansion  about  the  central  ray. 
Beatnwidth  in  this  method  is  the  distance  from  the  central 
ray  to  the  beam  Fresnel  volume  and  is  fully  determined  by 
equation  (59),  while  beamwidth  in  the  standard  Gaussian 
beam  method  is  usually  chosen  arbitrarily  and  not  given  any 
physical  meaning. 


3. 5. 4  Superposition  of  beams 

For  the  wavefield  obtained  by  superposition  of  all  beams  we 
shall  use  upper-case  U,  instead  of  lower-case  u,  which  is 
reserved  for  an  individual  beam.  Note  that  the  wavefield  u 
corresponding  to  an  individual  beam  is  a  function  of  vertical 
take-off  angle  and  azimuth  (6  and  <p),  which  specify  the 
central  ray  under  consideration.  Thus  we  shall  write 
u(S,  t,  6,  <p)  instead  of  u(S.  t ).  The  wavefield  U(S,  i)  will  be 
described  by  superposition  of  individual  beams  (rays)  with 
respect  to  6  and  p: 

U(S,  i)  =  f  f*u(S.  t,  6,  p)  d6  dp.  (66) 

When  the  integrand  of  equation  (66)  is  sufficiently  smooth 
for  a  given  S  and  /.  it  cun  be  discretized  as 

N  M 

ms,  <)=S  2  I,  t>r  <P*)A<5,A*-  (67) 

l-Ok-ll 

where  the  quantities  Ad,  and  Ad*  are  determined  from  a 
given  range  of  take-off  angles  6,  and  pk  (fervent  1983). 
The  wavefield  is  calculated  in  (66)  or  (67)  by  summing  up 
each  ray’s  contribution  at  a  specified  point,  its  wavelet 
having  a  Gaussian  distribution  both  in  space  and  time.  As  in 
the  Gaussian  beam  method  (Cerveny  1983).  this  method 
also  does  not  require  two-point  ray  tracing  to  compute  the 
seismic  wavefield. 

Since  the  energy  of  a  beam  is  conserved  along  the 
wavefront,  U(s,  l)  in  equation  (66)  can  be  rewritten  in  the 
wavefront  coordinates  by  using  the  Jacobians  in  equation 
(39) 


U(S- "  ■  LCdm”*  HW  -  (£r)1 

x  iA-(<J>y(7')(-/)*)  t/6  dp. 


(68) 


where  T  =  r  -  r.  As  a  shown  in  equation  (39).  the  Jacobians 
J i  and  J,  between  the  ray  centred  and  the  wavefront 
coordinates  give 

J,n,=n„  h,  =  R,a„  (69) 

where  or,  is  the  angular  difference  between  the  central  ray 
and  the  vicinity  ray  which  passds  through  at  point  (.s.  n ,)  in 
the  /-<?.  plane  (Fig.  3). 


3. 5. 5  Superposition  in  a  homogeneous  medium 

It  is  easy  to  demonstrate  that  the  superposition  integral 
returns  simple  ray  theory  in  a  homogeneous  medium.  In  a 
homogeneous  medium.  LHS,  t)  is  represented  simply 
because  Rt=  R:  =  S,  and  S  is  the  total  distance  from  the 
source  to  the  receiver.  The  parameters  in  a  homogeneous 
medium  are  given  as 

R,=S.  R, » ,V.  <!>=!, 


a  =  r,  +  A,  li  =  < ,  +  p. 

(70) 

and 

da  =  dd.  ilfi -dp. 

(71) 

where  c,  and  r.  are  constant 

Substituting  equations  (69), 
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(70).  (71)  into  (68).  then  gives 
\J2D  f* 


,  V2 D  r  f*  (Rtay  ,R2py  1 

1 0  V.tVF.F.  L J*, CXP  (j,  F, )  ~  (j3f)  J 


x  #,[y[T)(-i)k\da  dti 
V2  D 


(72) 


where  J  =J\J2  and  A:  =  0  in  a  homogeneous  medium, 
liquation  (72)  shows  that  the  displacement  of  U{S,  I )  in  a 
homogeneous  medium  is  proportional  to  l/.V,  the  distance 
between  source  and  receiver.  This  demonstrates  that  the 
superposition  is  a  high-frequency  solution  of  the  wave 
equation.  Note  that  the  superposition  is  independent  of  the 
choices  made  for  the  initial  conditions  for  the  spacing  of 
vicinity  rays. 


3.6  A  numerical  example  superposing  Gaussian  beams 
determined  fres  VRT 

A  numerical  examp'e  was  chosen  to  test  the  superposition  of 
Gaussian  beams  dehi.'d  from  the  VRT  system.  The  velocity 
model  (Fig.  S)  used  in  .he  test  is  the  same  as  that  used  in  the 
tests  of  the  accuracy  of  the  VRT  and  DRT  systems.  P 
velocity  is  continuous  h  this  model,  and  a  discontinuity  in  P 
velocity  gradient  exist,  at  10km  depth.  The  model  was  not 
designed  to  be  geophysically  realistic,  but  rather  to  illustrate 
the  theoretical  phenomena  near  a  caustic.  Fig.  10  shows  the 
results  of  ray  tracing,  indicating  a  triplication  in  the  range 
32.4-48.3  km  from  the  source.  Two  caustics  are  located  at 
*  =  32.4  and  46.3  km.  Fig.  11  shows  the  vertical  component 
synthetic  seismograms  computed  by  superposing  Gaussian 
beams  defined  from  vicinity  rays,  called  'VRT  seismograms'. 
To  calculate  these  seismograms,  64  rays  are  used  with 
Ad  =  P  successive  increments  in  take-off  angle.  An 
explosive  point-source  is  assumed,  and  the  initial  conditions 
are  ij,',  =  P  and  q„  =  0.  (Here,  the  choice  tj(,  =  Ad  is  made  to 
check  the  accuracy  of  the  VRT  computations  as  shown  in 


Figure  II.  The  vertical  component  Gaussian  beam  seismograms 
using  VRT  an  Fresnel  bcamwidths  for  the  model:  the  centre 
frequency  of  a  narrow-band  Gaussian  source  wavelet  is  5  Hz.  and 
the  receivers  are  located  at  the  surface  (z  =  ()  km). 

Figs  6  and  7.  To  improve  the  traveltime  estimation,  the 
choice  ?)o=  Ad/2  is  recommended.)  A  Gaussian  source 
(equation  55)  pulse  with  centre  frequency  5  Hz  is  used.  The 
beamwidths  are  taken  to  be  the  Fresnel  beamwidths  (59). 
WKBJ  synthetic  seismograms  are  computed  for  comparison 
with  the  VRT  seismograms  (Fig.  12).  As  shown  in  Figs  1 1 


Distance  (Kin) 

Figure  10.  Ray  trajectories  in  the  model  shown  in  Fig,  5:  the  triplication  zone  is  located  in  the  range  x  =  32.4-48.3  km. 
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Flgarc  12.  The  vertical  component  WKBJ  seismograms  for  the 
model  for  5  Hz:  the  conditions  are  the  same  as  in  the  VRT 
seismograms  in  Fig.  11  except  that  the  WKBJ  seismograms  were 
first  synthesized  for  a  delta  source-time  function  and  then  convolved 
with  a  narrow-hand  Gaussian  wavelet. 


and  12,  the  two  methods  closely  agree  with  one  another. 
Amplitude  differences  between  the  two  methods  are  less 
than  5  per  cent.  Diffractions  are  shown  near  the  caustics  in 
both  methods.  The  diffraction  near  the  caustic  at 
x  =  32.4  km  decays  faster  than  that  near  the  caustic  at 
x  «  46.4  km  because  the  beamwidth  (beam  Fresnel  volume) 
varies  more  slowly  at  the  former  than  at  the  latter  caustic. 
Some  differences  in  the  frequency  content  of  the  diffraction 
from  the  caustic  at  x»  32.4  km  can  be  seen.  These 
differences  were  generated  by  allowing  for  a  broad 
frequency  spectrum  in  the  WKBJ  method,  resulting  from  a 
delta  function  source  and  then  convolving  the  result  with  a 
narrow-band  source  pulse.  In  the  superposition  of  vicinity 
rays,  on  the  other  hand,  only  the  frequencies  contained  in 
the  narrow-band  source  pulse  are  considered,  because  the 
half-width  of  the  beam  Fresnel  volume  (beamwidth)  was 
computed  only  for  the  centre  frequency  of  the  narrow-band 
source  pulse. 

As  shown  in  Figs  II  and  12,  the  VRT,  the  WKBJ 
methods  produce  very  closely  matching  synthetic  seismo¬ 
grams.  The  minor  differences  among  the  methods  result 
from  the  characteristics  of  each  method.  Not  surprisingly, 
the  VRT  and  the  WKBJ  seismograms  are  nearly  alike 
because  both  methods  are  based  on  asymptotic  ray  theory. 


4  CONCLUSIONS 

The  VRT  system  traces  the  position  of  a  ray  in  the  ray 
centred  coordinate  system  of  a  reference  ray.  Although  the 
DRT  system  can  be  used  to  estimate  the  position  and 
traveltime  of  all  rays  in  the  paraxial  vicinity  of  the  reference 
ray,  it  breaks  down  rapidly  as  the  perpendicular  distance  q 
to  the  vicinity  ray  increases.  The  VRT  system  exhibits  a 
much  slower  breakdown  as  the  distance  to  the  vicinity  ray 
increases  and  closely  reproduces  the  ray  trajectory  obtained 
from  kinematic  ray  tracing  even  at  large  values  of  q. 
Quantities  obtained  from  integrating  the  VRT  system  can  be 
used  to  determine  geometric  spreading,  wavefront  curva¬ 
ture,  and  KMAH  index  with  greater  numerical  precision 
than  is  possible  by  differencing  kinematically  traced  rays. 

The  VRT  system  does  not  require  transformation  from 
ray  cenvrcd  to  ray  coordinates.  A  reduction  in  the  number 
of  required  coordinate  systems  in  VRT  leads  to  a  fewer 
number  of  equations  needed  to  specify  quantities  related  to 
the  wavefront.  The  VRT  system  is  specified  by  only  four 
equations  in  addition  to  the  kinematic  ray  tracing  equations. 
By  contrast,  the  standard  dynamic  ray  tracing  equations 
based  on  paraxial  approximations  require  eight  equations. 
Unlike  the  standard  DRT  system,  the  VRT  system  does  not 
require  the  evaluation  of  second  spatial  derivatives  of 
velocity  along  a  ray.  The  VRT  equations  will  thus  also  have 
advantages  over  standard  DRT  in  simplifying  the  para- 
metrization  of  the  medium.  Since  only  first  spatial 
derivatives  of  velocity  are  used  in  the  VRT  system,  no  phase 
matching  is  required  at  discontinuities  in  velocity  gradient. 
At  velocity  discontinuities,  new  initial  conditions  on  the 
vicinity  rays  are  determined  by  applying  Snell’s  law  to  the 
central  ray  and  the  vicinity  rays  separately. 

The  quantities  obtained  from  VRT  may  be  easily 
converted  to  those  obtained  from  DRT  and  vice  versa.  Since 
VRT  provides  a  much  more  accurate  prediction  for  the 
trajectory  of  a  ray  in  the  vicinity  of  a  reference  ray,  the 
relations  between  VRT  and  DRT  can  be  used  to  estimate 
the  'paraxial  vicinity'  in  which  errors  associated  with  the 
paraxial  approximations  of  DRT  are  small. 

In  either  DRT  or  VRT  a  physical  definition  of  beamwidth 
may  be  obtained  from  the  concept  of  a  beam  Fresnel 
volume.  An  example  calculation  demonstrates  that  this 
definition  of  beamwidth  can  approximate  diffraction  effects. 
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ABSTRACT 

Vicinity  ray  tracing  (Section  I)  is  used  to  determine  the  focusing/defocusing  and  multipathing  in¬ 
duced  by  the  Aleutian  slab  on  the  P  waves  radiated  by  Aleutian  nuclear  tests.  The  geographic 
pattern  of  amplitude  anomalies  predicted  from  slab  models  proposed  from  modeling  of  P  travel 
time  anomalies  is  consistent  with  the  observed  geographic  distribution  of  m*,  anomalies  from  the 
Aleutian  tests  and  the  results  of  studies  of  P  waveforms  observed  from  shallow  focus  earthquakes 
occuring  within  the  Aleutian  island  ridge.  A  broad  zone  of  reduced  P  amplitudes  is  predicted  at 
northerly  azimuths  at  distances  greater  than  70°.  This  shadow  zone  is  also  likely  to  be  associated 
with  pulse  broadening  of  long  period  and  broadband  waveforms.  A  network  average  of  mt,  for 
shallow  focus  events  on  Amchitka  Island  is  predicted  to  underestimate  the  true  size  of  the  event  by 
as  much  as  0.4  mi,  units  if  all  stations  in  the  network  average  are  located  within  the  shadow  zone  of 
the  slab.  For  the  most  probable  distribution  of  stations  used  in  mj,  estimates  of  the  Aleutian  tests, 
this  bias  is  reduced  to  a  negative  bias  of  0.1  mj,  units.  Due  to  the  large  areal  overlap  of  European 
and  Canadian  stations  with  the  slab  shadow  zone,  weighting  of  stations  by  focal  sphere  solid  angle 
in  a  network  average  cannot  reduce  this  negative  bias  much  below  0.07  vii,  units.  The  documented 
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existence  of  a  similar  defocusing  anomaly  beneath  a  portion  of  NTS,  however,  may  eliminate  the 
need  of  making  a  correction  for  relative  focusing  and  defocusing  when  comparing  mj’s  of  Aleutian 
tests  with  those  of  NTS  tests. 

INTRODUCTION 

The  object  of  this  study  is  to  determine  the  yield  bias  of  underground  nuclear  tests  induced  by 
the  presence  of  a  high  velocity  descending  slab  beneath  the  test  site.  Specifically,  investigation  is 
made  of  the  effect  of  the  Aleutian  slab  on  the  measured  mj  of  the  underground  tests  Longshot,  Milro, 
and  Cannikan  P  wave  seismograms  are  synthesized  using  dynamic  ray  tracing  and  superposition 
of  Gaussian  beams  in  three-dimensional  models  of  the  Aleutian  slab  determined  from  P  travel 
time  delays.  Focusing  and  defocusing  and  multipathing  at  teleseismic  distances  are  evaluated  by 
comparison  of  observed  with  synthetic  seismograms  of  the  Aleutian  tests. 

COMPUTATIONAL  METHODS  AND  SLAB  MODELS 

Theoretical  amplitudes  and  travel  times  were  computed  using  vicinity  ray  tracing  (Section  I) 
in  the  long  slab  model  of  Boyd  and  Creager  Crcager  (1991)  for  the  Aleutian  slab  using  source 
positions  corresponding  to  underground  nuclear  tests  in  the  island  ridge  adjacent  to  the  slab.  In 
constructing  the  model  shown  in  figure  1,  the  raw  thermal  model  was  obtained  from  Crcager  (per¬ 
sonal  communication)  and  converted  to  a  P  velocity  model  by  assuming  the  temperature  derivative 
of  P  velocity  used  by  Boyd  and  Creager,  dVp/dT  =  0.5ms-1  A’-1.  Details  of  the  amplitude  calcu¬ 
lation  are  described  in  Scientific  Report  No.  2  of  this  project,  which  also  discusses  the  results  of 
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experiments  fore  several  different  types  of  slab  models  and  variations  of  source  position  relative  to 
the  slab. 


RESULTS 

Figure  2  shows  P  and  PcP  rays  predicted  for  Amchitka  tests  in  the  Boyd  and  Creager  model. 
The  rays  are  shown  for  a  2-D  cross  section,  perpendicular  to  the  strike  of  the  '•’ab.  Multipathing 
can  be  observed  at  the  great  circle  distances  42°  to  53°  for  P  waves  and  around  12°  for  PcP  waves. 

Amplitude  and  travel  times  were  calculated  in  models  with  and  without  the  slab  using  PREM 
as  a  reference  model.  P  amplitude  anomalies  are  shown  in  figure  3,  contoured  in  rrn  residuals. 
A  geographic  plotting  convention  is  used  rather  than  a  focal  sphere  plot.  The  epicenter  is  at  the 
center  of  the  sphere,  the  inner  circle  corresponds  to  the  area  at  distances  less  than  or  equal  to  35°, 
the  outer  circle  corresponds  to  core  grazing  distances. 

Highest  amplitudes  are  predicted  to  occur  around  an  annular  region  at  42°  to  53°.  Low  ampli¬ 
tudes  are  predicted  to  occur  at  longer  distances  outside  this  ring,  nearly  everywhere  on  the  dipping 
side  of  the  slab.  Peak  mj  residuals  in  the  observed  data  (figure  4)  are  bounded  by  0.6  m&  units.  The 
mb  anomalies  in  the  predictions  (figure  3)  are  uniformly  positive  (amplitudes  are  focused  relative 
to  those  observed  in  a  reference  model).  If  a  baseline  of  0.35  m&  is  subtracted  from  the  predictions 
shown  in  Figure  3,  then  the  strength  of  the  predicted  anomaly  pattern  is  smaller  than  the  anomaly 
pattern  shown  in  the  data  by  nearly  a  factor  of  2. 

The  predicted  variation  of  up  to  0.35  in  logioA,  however,  does  agree  with  the  predictions  of 
Sleep  (1973),  who  predicted  a  factor  of  2.5  variation  in  teleseismic  P  amplitudes  of  shallow  focus 
events  in  the  Aleutian  ridge.  Sleep  corrected  for  receiver  effects  by  dividing  measured  amplitudes 
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of  Aleutian  ridge  events  by  amplitudes  recorded  by  the  same  stations  for  nearby  events  that  are 
largely  unaffected  by  slab  focusing  and  defocusing.  The  higher  intensity  of  ms  variations  shown  in 
Figure  4,  consistent  with  nearly  a  factor  of  10  in  amplitude,  likely  contains  a  significant  component 
of  variation  due  to  the  effects  of  receiver  structure  and  regional  variation  of  attenuation.  The  much 
closer  agreement  with  amplitude  anomalies  corrected  for  these  effects  suggests  that  most  raw  ms 
data  will  have  large  effects  due  to  both  receiver  structure  and  regional  variations  in  path  attenuation 
as  well  near  source  focusing  and  defocusing.  The  unmodelled  receiver  and  path  effects  present  in 
the  ms  data  shown  in  figure  4  are  not  large  enough  to  obscure  a  slab  shadow  zone  visible  in  the 
concentration  of  negative  mb  anomalies  at  northerly  azimuths. 

The  results  shown  here  for  a  long  slab  model  are  essentially  the  same  as  those  discussed  in  the 
studies  of  Sleep  (1973)  and  Davies  and  Julian  (1972)  for  slab  models  extending  to  shallower  depths. 
Predicted  ms  anomalies  for  Aleutian  ridge  events  are  thus  not  sensitive  to  details  of  slab  structure 
below  650  km  depth. 

If  stations  used  in  a  network  averaged  mb  lie  entirely  within  the  slab  shadow  zone,  the  estimated 
mb  would  be  up  to  0.4  ms  units  smaller  than  that  estimated  if  the  slab  structure  were  not  present. 
This  bias  drops  to  a  -0.1  ms  if  the  network  shown  in  figure  4  is  used  to  construct  ms.  The 
values  plotted  in  figure  4  are  taken  from  studies  by  McLaughlin  et  al.  (1990),  in  which  clipped 
and  poor  signal  to  noise  waveforms  used  in  a  maximum  liklihood  estimate  are  removed.  The 
regional  concentrations  of  stations  is  representative  of  the  most  probable  network  used  in  a  network 
average  of  Aleutian  events,  large  gaps  in  coverage  corresponding  to  oceans  and  poorly  instrumented 
continents.  Weighting  of  stations  by  focal  angle  cannot  reduce  the  intensity  of  the  magnitude  bias 


much  below  0.1  mj  units  because  much  of  this  network  overlaps  a  large  areal  extent  of  the  slab 
shadow  zone. 

Calculations  of  body  wave  amplitudes  and  travel  times  were  also  made  using  earthquake  sources 
in  the  Aleutian  slab  and  slabs  in  other  regions  (Scientific  Report  No.  2  of  this  project).  For  earth¬ 
quakes  located  within  slabs,  it  was  assumed  that  regions  of  defocusing  correspond  to  regions  of 
maximum  broadening  and  complexity  in  body  waves.  Slabs  that  thicken  or  have  a  reduced  velocity 
contrast  below  650  km  depth  predict  a  different  regional  pattern  of  waveform  broadening  compared 
to  that  predicted  by  slabs  that  penetrate  the  650  km  discontinuity  for  a  long  distance  as  a  thin 
tabular  structure.  Data  from  the  Kuril-Kamchatka  slab  are  consistent  with  advective  thickening 
or  reduced  velocity  contrast  below  650  km  depth.  The  particular  pattern  of  S  and  ScS  waveform 
broadening  in  North  America  from  deep  focus  events  in  this  region  is  more  likely  to  be  a  conse¬ 
quence  of  a  slab  effect  than  an  attenuation  effect. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Slab  Effects  on  mi,  of  Aleutian  Tests 


The  sparseness  and  noise  in  raw  data  (figure  4)  do  not  allow  any  more  definitive  conclusion 
to  be  made  than  the  concentration  of  negative  mj,  anomalies  to  the  north  at  distances  greater  than 
53°  are  generally  consistent  with  a  shadow  zone  predicted  both  in  the  ray  tracing  shown  in  figure 
2  and  the  synthetic  anomaly  pattern  in  figure  3.  If  mi,  data  are  corrected  for  receiver  and  path 
attenuation  effects,  however,  the  agreement  between  observed  and  predicted  log  amplitude  data 
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is  excellent  for  the  focusing  and  defocusing  effects  of  the  Aleutian  slab.  A  negative  bias  of  0.1 
mb  units  is  predicted  from  the  most  probable  network  average  of  m&.  A  similar  defocusing  body 
lies  beneath  the  Pahute  Mesa  region  of  NTS  (Taylor,  1983),  having  a  similar  relative  location  of 
teleseismic  rays  reaching  a  dense  concentration  of  stations  in  a  global  network  (Cormier,  1987). 
This  situation  would  will  tend  to  mitigate  the  need  for  or  reduce  the  value  of  a  relative  correction 
for  near  source  focusing  and  defocusing  when  NTS  and  Aleutian  mj’ s  are  compared. 

The  lowest  amplitudes  in  predicted  and  observed  data  are  due  north,  perpendicular  to  the 
strike  of  the  Aleutian  slab,  at  ranges  exceeding  53°.  This  is  the  region  in  which  evidence  of 
pulse  broadening  has  been  reported  in  long  period  and  broadband  waveforms  from  shallow  focus 
Aleutian  events  (Engdahl  et  al.,  1989).  The  coincidence  of  the  region  of  waveform  broadening  with 
the  region  of  defocusing  suggests  that  the  pulse  broadening  is  caused  by  slab  diffraction  (Vidalc, 
1987;  Cormier,  1989). 

Methods  of  Waveform  Synthesis  for  Slab  Effects 

New  constraints  can  be  derived  from  digital  waveform  data  through  investigations  of  the  effect 
of  the  Aleutian  slab  on  P  waveform  broadening  and  complexity.  Both  finite  difference  (Vidale, 
1987)  and  Gaussian  beam  synthetics  (Cormier,  1989;  Weber,  1990)  have  been  applied  to  calculate 
these  efiV’ts.  The  computational  expense  of  the  finite  difference  approach  practically  limits  its  use 
to  two-dimensional  geometry.  Although  limited  three-dimensional  results  have  been  obtained  with 
the  standard  Gaussian  beam  approach,  it  is  close  to  the  limits  of  its  applicability  with  most  slab 
models  and  cannot  properly  describe  the  Fresnel  zone  in  the  vicinity  of  the  source  region  in  the  slab 


(Cormier  and  Kim,  1990).  A  promising  new  approach  to  calculating  frequency  dependent  effects  of 
slab  structure  is  a  form  of  Born  scattering  approximation  that  does  not  require  a  reference  medium 
(Coates  and  Chapman,  1991).  This  approach  should  be  capable  of  including  the  low  frequency 
scattering  of  energy  into  the  slab  diffracted  phase  by  the  zones  of  high  gradient  defining  the  slab. 
Contours  of  Fresnel  zones  can  be  calculated  from  intermediate  results  of  the  Born  scattering  theory. 
Figure  6  shows  the  results  of  a  Fresnel  zone  calculation  for  source  within  the  high  velocity  slab  shown 
in  figure  5.  Combined  with  maps  of  zones  of  strong  gradient  (figure  7),  the  Fresnel  zone  maps  can 
be  used  to  qualitatively  predict  which  source/receiver  pairs  will  exhibit  waveform  broadening. 
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Figure  1:  P  velocity  model  of  the  Aleutian  slab  from  Creager  and  Boyd  (1989).  Location  of 
Amchitka  nuclear  tests  relative  to  the  slab  structure  is  shown  by  the  asterisk  at  the  top. 


Figure  2:  P  and  PcP  ray  paths  for  the  slab  model  and  source  position  shwon  in  Figure  1.  Cross 
section  of  the  Earth  is  in  a  plane  perpendicular  to  the  strike  of  the  slab,  containing  the  source. 


Figure  3:  Predicted  P  amplitude  anomalies  contoured  as  magnitude  residuals  for  the  slab  and  source 
position  shown  in  Figure  1.  Plot  is  an  equal  area  geographic  plot,  with  North  at  the  top.  The 
blank  region  in  the  center  corresponds  to  distances  less  than  or  equal  to  35°  great  circle  degrees. 
The  outer  radius  corresponds  to  the  distance  of  core  grazing  P  waves. 
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Figure  4:  Observed  mb  residuals  of  the  Cannikin  nuclear  test  calculated  with  respect  to  the  maxi¬ 
mum  likelihood  mj,  reported  by  McLaughlin  (1990) 
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Figure  5:  Shear  velocity  contours  of  a  deeply  penetrating  slab  model  having  advective  thickening 
below  650  km. 
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arival  time  of  scattered  low  frequency  energy  referenced  to  the 
deep  focus  source  located  in  the  slab  model  shown  in  Figure  5. 
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Figure  7:  Contours  of  | W\/V  for  shear  velocity  of  the  slab  model  shown  in  Figure  5.  The  origin 
of  scattered  low  frequency  energy  in  a  slab  diffracted  pulse  can  be  found  by  overlaying  Figure  6 
onto  this  figure. 
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ABSTRACT 

Virtually  all  regional  phases  can  be  strongly  affected  by  vertical  velocity  gradients.  The  best 
known  effects  are  on  the  Pn  and  Sn,  in  which  small  changes  in  the  vertical  velocity  gradient 
beneath  the  Moho  produce  large  changes  in  the  decay  of  Pn  and  Sn  with  distance.  Methods  of 
synthesizing  complete  regional  seismograms  often  inadvertently  ignore  the  effect  of  crustal  gradients 
by  parameterizing  the  Earth  model  with  thick,  planar  homogeneous  layers.  To  address  this  problem 
we  have  modified  the  locked  mode  method  of  synthesizing  complete  regional  seismograms  to  include 


the  Langer  uniform  asymptotic  approximation  to  vertical  wavefunctions  within  layers  having  linear 
vertical  velocity  gradients.  Synthesis  of  complete  regional  seismograms  using  the  Langer-locked 
mode  confirm  that  the  Pn  and  Sn  phases  are  strongly  affected  by  the  magnitude  of  the  velocity 
gradients  beneath  the  Moho,  but  that  Lg  is  only  weakly  affected  by  the  details  of  crustal  layering. 

Tests  were  made  to  quantify  the  error  in  the  use  of  the  Langer  approximation  as  the  magnitude 
of  the  vertical  gradient  increases  and/or  frequency  decreases.  At  sufficiently  small  magnitude  of 
gradient  and/or  high  frequency,  good  agreement  can  be  obtained  between  synthetics  computed 
using  the  Langer-locked  mode  method,  the  colocation  method,  and  the  conventional  locked  mode 
method  in  models  parameterized  by  thin  homogeneous  layers.  Errors  introduced  by  the  use  of  the 
Langer  approximation  in  calculated  pole  positions,  residues,  and  eigenfunctions  are  bounded  by  5% 
for  frequencies  /  >  5  |VV|.  An  upper  bound  to  the  error  in  the  time  domain  can  be  estimated  from 
this  inequality  using  either  the  peak  frequency  in  a  narrow  pass  band  or  the  lowest  frequency  in  a 
broad  pass  band.  When  10  or  more  thin  homogeneous  i.  yers  are  required  to  represent  accurately 
the  seismic  wavefield  in  a  gradient  layer,  it  is  usually  more  efficient  to  represent  the  gradient  layer 
by  continuously  varying  functions  in  the  vertical  direction  and  employ  the  Langer  approximation, 
provided  the  errors  in  the  Langer  approximation  remain  within  acceptable  limits.  By  reducing  the 
number  of  parameters  needed  to  describe  an  earth  model,  the  Langer-locked  mode  method  simplifies 
the  inverse  problem  of  determining  structure  using  observed  and  synthetic  complete  seismograms. 
It  also  facilitates  the  use  of  known  relations  for  the  effects  of  continuously  varying  pressure  and 
temperature  on  elastic  moduli  and  density. 
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INTRODUCTION 


Complete  seismograms  at  local  and  regional  distances  are  now  routinely  computed  in  plane  layered 
models  for  a  variety  of  source  receiver  geometries,  source  depths,  and  source  types  by  integrating  or 
summing  over  wavenumbers  (e.g.  Bouchon  and  Aki,  1977;  Kind  1978;  Wang  and  Hermann,  1980; 
Mandal  and  Mitchell,  1986;  Mandal  and  Toksoz,  1990)  or  summing  locked  or  leaky  modes  (e.g. 
Harvey,  1981;  Kerry,  1981;  Haddon,  1986;  Nolct  et  al.,  1989).  The  computational  expense  of  these 
calculations  remains  relatively  cheap  as  long  as  the  crust  and  upper  mantle  model  can  be  described 
by  a  small  number  of  planar,  homogeneous  layers. 

Seismograms  synthesized  in  models  composed  of  a  small  number  of  plane  homogeneous  layers 
ignore  the  continuous  depth  dependence  of  clastic  moduli.  Usually  seismograms  are  synthesized 
in  simple  models  composed  of  two  or  three  homogeneous  layers  of  crust  overlying  a  homogeneous 
lid,  low  velocity  zone,  and  upper  mantle  beneath  the  lid.  Since  Earth  curvature  is  ignored  in  these 
calculations,  the  model  is  effectively  one  in  which  each  layer  has  a  small  negative  gradient  with 
depth. 

A  truly  realistic  earth  model  would  include  three-dimensional  distributions  of  heterogeneities 
having  a  broad  spectrum  of  scale  lengths.  At  wavelengths  much  larger  than  the  largest  scale  length 
of  heterogeneity,  the  propagation  of  the  wavefield  in  such  a  model  can  be  accurately  calculated  in 
an  equivalent  model,  which  is  generally  anisotropic.  If  there  are  no  preferred  orientations,  shapes, 
or  periodicities  in  the  distribution  of  heterogeneities,  and  if  the  ensemble  averages  of  densities  and 
elastic  moduli  primarily  increase  or  decrease  as  a  function  of  depth,  then  the  equivalent  model  can  be 
most  simply  represented  by  a  small  number  of  vertically  inhomogeneous,  isotropic  layers.  Synthesis 
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of  the  wavefield  in  this  long  wavelength,  equivalent  model  ignores  the  effects  of  scattering,  which 
become  important  as  the  wavelength  approaches  the  size  of  the  scale  length  of  heterogeneities, 
but  still  includes  important  effects  of  mean  vertical  gradients  on  the  interference  seismic  phases 
propagating  horizontally  between  vertical  discontinuities.  Such  a  model  also  facilitates  comparison 
of  estimated  velocities  and  densities  with  the  predictions  of  theoretical  and  empirical  relations  for 
the  effects  of  vertically  varying  pressure,  temperature,  pore  fluids,  and  crack  density. 

Virtually  all  of  the  regional  phases  can  be  strongly  affected  by  vertical  velocity  gradients.  One 
example  is  dispersion  of  the  fundamental  mode  Rayleigh  wave,  or  Rg  phase,  at  local  and  regional 
distances.  Although  often  strongly  scattered  by  surface  topography  and  near  surface  heterogeneity, 
good  examples  of  dispersed  Rg  wavetrains  can  sometimes  be  found  at  epicentral  distances  on  the 
order  of  50  km  (e.g.,  Kafka  and  Reiter,  1987).  Simple  crustal  models  having  a  homogeneous  layer 
of  10  km  or  more  thickness  at  the  surface  produce  an  unrealistically  impulsive,  undispersed  Rg 
arrival.  To  reproduce  the  observed  Rg  dispersion  properly,  the  velocity  and  density  model  must 
have  a  strong  positive  gradient  with  increasing  depth. 

Perhaps  the  best  known  effects  of  gradients  on  regional  phases  are  those  on  the  Pn  and  Sn 
phases.  In  a  plane  layered  model  composed  of  a  crust  overlying  a  thick  upper  mantle  lid,  the  Pn 
and  Sn  phases  are  .  'assical  hcadwaves  traveling  just  beneath  the  Moho.  Rill  (1971)  and  Cerveny 
and  Ravindra  (1971)  have  shown  how  gradients  transform  classical  headwaves  into  interference 
headwaves  or  ’’whispering  gallery  waves”  (e.g.,  Cormier  and  Richards,  1977;  Menke  and  Richards, 
1980).  The  distance  decay  of  both  classical  and  interference  headwaves  is  frequency  dependent. 
The  interference  headwave  decays  much  more  slowly  with  distance  than  the  classical  headwave 
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for  a  mantle  having  a  positive  gradient  with  depth  below  the  Moho  and/or  for  a  Moho  boundary 
incorporating  earth  curvature. 

In  this  paper,  we  describe  experiments  in  which  vertical  velocity  gradients  are  incorporated  in 
models  of  the  crust  and  upper  mantle  using  the  locked  mode  method  together  with  a  high  frequency, 
asymptotic  approximation  to  the  vertical  wavefunctions  (Appendix  I)  in  each  vertically  inhomo¬ 
geneous  layer.  The  purpose  of  these  experiments  is  to  (1)  quantify  the  breakdown  in  asymptotic 
(Langer)  approximation  to  the  vertical  wavefunctions  as  frequency  decreases  and/or  the  magnitude 
of  the  vertical  gradient  increases  and  to  (2)  illustrate  and  review  some  of  the  effects  of  vertical 
velocity  gradients  on  the  propagation  of  regional  seismic  phases.  Although  all  of  the  example 
seismograms  are  synthesized  for  high  frequency,  local  and  regional  seismic  phases,  the  programs 
and  algorithms  may  also  be  applied  to  lower  and  intermediate  frequency  phases  at  longer  ranges, 
overlapping  the  high  end  of  the  frequency  band  traditionally  included  in  normal  mode  studies  of 
the  whole  earth. 

The  paper  begins  with  a  brief  review  of  the  locked  mode  method.  Mathematical  details  of  the 
Langer  approximation  and  its  incorporation  in  the  locked  mode  method  are  described  in  Appendices 
I  and  II.  The  body  of  the  paper  describes  the  results  of  tests  conducted  to  determine  the  accuracy 
of  the  Langer  approximation  and  how  it  breaks  down  as  the  gradient  in  the  layer  increases.  A 
discussion  and  example  show  how  depth  and  frequency  dependent  attenuation  can  be  included  in 
the  Langer-locked  mode  method.  The  paper  concludes  with  examples  of  synthetic  seismograms  in 
the  100  to  1000  km  range,  showing  that  gradients  near  the  free  surface  and  Moho  can  radically 
affect  the  propagation  of  some  of  the  principal  regional  phases. 


37 


REVIEW  OF  THE  LOCKED  MODE  METHOD 


Representation 

Following  Harvey  (1981;  1985),  the  complex  displacement  spectra  are  evaluated  from 
Ru{u,xn9T,zr)  =  -RcI-Rfil-'Ym  «A(n>w)  RE(n,u,za)  ^(n)w,m,ir)Or,:r) 


Lu(u),xr,0rjZr)  =  ~LfiI  -  £,A(n,«)  LST(n,u,m)  LE(n,w,z.)  L^(n,w,m,xr,^,zr) 

n  m 

where  the  subscripts  R  and  L  denote  Rayleigh  and  Love  modes  respectively. 
flA  and  iS.  area  scalar  amplitude  factors  defined  by 


flA(n,w)  =  _ 


RkY,  3(0) 

dYn(0)/dk 


£,A(n,w)  =  - 


LkP2(0) 

dDi{0)/dk 


lij(O)  and  2?,(0)  are  evaluated  at  the  free  surface  (z  =  0 )  and  are  defined  in  Appendices  I  and 
II.  r2t  and  £,Sr  are  row  vectors  determined  from  the  source  jump  vectors.  Rtp  and  ixj)  are  defined 
from  products  of  eigenfunctions  for  displacement  (rEi,  rEi)  and  iE\  evaluated  at  the  receiver 
depth  ( z  =  zr )  and  vector  cylindrical  harmonics  P,  B,  C: 

R$(n,u>,m,Xr,Or,Zr)  =  rEi (n, u>,  *r)P(n,  tv,  m,  xr,  0r)  +  R^n,  w,  *r)6(n,  W,  m,  arr,  0r) 


Llp(n,U,m,Xr,Or,Zr)  =  /,£?!  (n,  U>,  Zr)C(7l,  W,  m,  Xr,  dT) 
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Ha  /  and  i0I  arc  branch  cut  integrals,  which  account  for  energy  that  cannot  be  represented  by 
normal  modes,  and  are  associated  with  near  vertically  propagating  P  and  S  waves  that  leak  into 
the  halfspace.  The  locked  mode  method  does  not  evaluate  the  branch  cut  integrals.  It  chooses 
the  halfspace  to  be  sufficiently  deep  and  fast  such  that  all  of  the  energy  important  to  a  particular 
time  window  at  a  particular  distance  can  be  accurately  represented  by  the  locked  mode  summation 
alone. 

Implementation  with  the  Langer  Approximation 

Seismograms  arc  synthesized  by  evaluating  the  complex  spectra  at  discrete  frequencies  and  invert¬ 
ing  to  the  time  domain  by  a  fast  Fourier  transform.  A  small  complex  frequency  can  be  added  to 
attenuate  all  arrivals  outside  of  the  finite  time  window  given  by  the  folding  frequency  of  the  discrete 
Fourier  transform  (Rosenbaum,  1974;  Muller  and  Schott,  1981).  Harvey  (1981, 1985)  gives  detailed 
derivations  of  the  locked  mode  method  and  describes  its  implementation  in  media  described  by 
homogeneous  layers.  The  principal  modifications  of  the  method  for  use  with  the  Langer  approxi¬ 
mation  are  concerned  with  the  calculation  of  the  eigenfunction  vector  E  and  the  scalar  amplitude 
factors  a  A  and  £,A  (Appendices  I  and  II)  The  partial  derivatives  with  respect  to  k  appearing  in  the 
amplitude  factors  are  calculated  by  difference  derivatives. 

Since  the  Langer  approximation  describes  turning  rays  in  each  layer,  it  is  desirable  to  include 
the  effects  of  earth  curvature  using  a  locked  mode  representation  in  fully  spherical  rather  than 
cylindrical  coordinates.  Accordingly,  all  formulae  in  Appendices  I  and  II  are  given  as  functions  of 
radially  varying  velocities  and  density  and  ray  parameter,  p,  in  a  spherical  earth.  For  applications 
to  high  frequency  regional  seismograms,  insignificant  differences  will  be  obtained  between  results 
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using  a  representation  written  in  terms  of  cylindrical  harmonics  P,B,C  versus  the  results  using  a 
representation  written  in  terms  of  spherical  harmonics,  provided  that  the  wavenumber  k  is  associ¬ 
ated  with  the  ray  parameter  p  in  a  spherical  earth  by  the  relation,  k  =  u>p/Re ,  where  Rt  is  the 
radius  of  the  earth.  For  values  of  the  non-dimensional  products  kx  or  up  exceeding  5,  the  correction 
factor  needed  to  make  cylindrical  harmonics  reproduce  the  result  obtained  with  spherical  harmonics 
is  given  by  A/ sin(A)  (Muller,  1971;  Cormier  and  Richards,  1989).  Since  this  factor  varies  only 
between  1.0  and  1.002  for  distances  between  0  and  1000  km.,  we  have  omitted  this  correction  in 
our  test  calculations  and  simply  used  the  cylindrical  harmonic  representation  of  equations  1  and 
3. 

The  Langer  approximation  can  also  be  implemented  in  methods  of  synthesizing  complete  seis¬ 
mograms  that  numerically  integrate  over  horizontal  wavenumber  and  slowness  (Cormier,  1980). 
The  primary  advantage  of  the  locked  mode  method  is  that  most  of  the  computational  effort  in¬ 
volved  in  the  calculation  of  the  amplitude  factors  and  eigenfunctions  can  be  cataloged  for  use  with 
different  source-receiver  geometries  and  different  moment  tensor  representations  of  point  sources. 
Although  response  functions  can  be  similarly  cataloged  in  approaches  that  integrate  or  sum  over 
wavenumber  or  slowness,  this  is  rarely  done  in  practice.  A  secondary  advantage  of  the  locked  mode 
method  is  that  a  large  body  of  literature  exists  in  modal  notation  on  inversion  for  structure  and 
source  parameters.  The  analysis  of  problems  using  normal  modes  of  the  whole  Earth  at  low  fre¬ 
quency  and  long  range  can  usually  be  directly  adapted  to  higher  frequency  and  shorter  range  using 
locked  modes  (e.g.,  Gomberg  and  Masters,  1988). 
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ACCURACY  OF  THE  LANGER  APPROXIMATION 


•  > 

H 

/ 

The  Scalelength  Parameter  and  Reference  Models 

The  Langer  approximation  assumes  decoupling  between  P  and  S  waves  and  up-  and  down-going 
waves  in  each  gradient  layer,  and  the  criteria  for  its  accuracy  are  thus  similar  to  those  used  in  ray- 
asymptotic  solutions  to  the  elastodynamic  equation  of  motion  in  inhomogeneous  media  (Richards, 
197G).  Qualitatively,  the  Langer  as  well  as  all  other  ray-asymptotic  approximations  to  the  solution 
of  the  elastodynamic  equation  of  motion  are  known  to  become  less  accurate  as  non-dimensional 
ratios  A/(v/|Vu|)  increase,  where  v  is  a  P  or  S  velocity  or  density  (Richards,  1976;  Chapman,  1974). 
Another  way  in  which  this  is  commonly  phrased  is  that  the  wavelength  must  be  much  smaller 
than  the  scalelength  of  the  medium,  L ,  where  L  is  the  maximum  of  (a/|Va|,/?/|V/?|,/>/|Vp|) 
(Beydoun  and  Ben-Menahem,  1985).  A  goal  in  this  study  is  to  quantify  the  breakdown  in  the 
Langer  approximation  as  the  scalelength  of  gradient  layers  decrease,  determining  exactly  how  large 
the  ratio  X/L  can  be  before  errors  in  calculated  displacement  exceed  some  specified  bound. 

The  first  step  in  such  a  study  is  to  choose  accurate  reference  synthetic  seismograms  in  models 
having  strong  gradients.  Spudich  and  Ascher  (1983)  published  synthetic  seismograms  calculated 
by  the  numerical  colocation  method  for  a  simple  model  consisting  of  a  gradient  over  halfspace. 
The  gradient  layer  in  this  model  was  parameterized  by  a  sequence  of  40  thin  layers  (Figure  1),  the 
width  of  each  thin  layer  approximately  equal  to  one- tenth  the  wavelength  of  shear  waves  at  1  Hz. 
Excellent  agreement  was  found  between  the  locked  mode  synthetics  and  the  colocation  synthetics. 
This  result  confirmed  that  locked  mode  synthetics  computed  in  models  in  which  gradient  layers 
are  represented  by  thin  layers  can  be  used  as  accurate  reference  synthetics  to  test  the  Langer 
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approximation. 

To  test  the  accuracy  of  the  Langer  approximation,  seismograms  were  synthesized  using  the 
locked  mode  method  using  the  Langer  approximation  in  a  series  of  models  with  increasing  gradients 
in  P  and  S  velocity  and  density  in  a  layer  over  a  halfspace  (Figure  2).  The  series  of  models  and 
the  frequency  band  of  synthesis  were  purposely  chosen  to  explore  the  breakdown  of  the  Langer 
approximation  at  large  of  X/L. 

Dispersion  Curves 

Figure  3  compares  the  dispersion  curves  of  the  locked  Love  and  Rayleigh  modes  calculated  with 
Langer  approximation  in  a  thick  continuous  gradient  layer  ( dashed  line)  with  those  calculated  by 
parameterizing  the  gradient  layer  with  thin  homogeneous  layers  ( solid  line).  Results  are  shown  only 
for  model  1,  which  has  the  strongest  surface  gradient  in  Figure  2.  Even  for  the  severe  gradients  in 
the  surface  layer  of  model  1,  in  which  X/L  ranges  from  40  at  0.1  Hz  to  3  at  1  Hz,  dispersion  curves 
calculated  using  the  Langer  approximation  closely  track  those  calculated  in  the  model  parameterized 
by  thin  homogeneous  layers.  As  expected,  the  errors  in  the  Langer  approximation  are  generally 
largest  at  lower  frequency,  where  X/L  is  largest.  The  primary  region  of  error  occurs  for  the  low 
frequencies  of  the  fundamental  mode.  This  is  not  unexpected  since  most  of  the  energy  of  the 

fundamental  mode  in  this  frequency  band  is  confined  to  the  strong  gradient  layer  near  the  surface. 
Errors  in  the  dispersion  of  Love  modes  show  a  simple  dependence  on  frequency  and  begin  to  be 
acceptably  small  enough  for  accurate  calculations  even  at  values  of  X/L  approaching  1.  Errors 
in  the  dispersion  of  Rayleigh  modes  have  a  slightly  more  complex  behavior.  In  the  dispersion  of 
Rayleigh  modes,  note  a  region  of  phase  velocities  and  frequencies  between  4.2  to  6.0  km/sec  and 
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0.4  to  0.7  Hz.  In  this  region  of  the  dispersion  plot,  the  thin  layered  calculation  shows  a  group  of 
4  to  6  dispersion  curves  having  two  inflections  in  curvature  occurring  over  short  ranges  of  about 
0.2  km/sec  in  velocity  and  0.1  Hz  in  frequency.  The  Rayleigh  dispersion  curves  calculated  by 
the  Langer  approximation  do  not  have  these  inflection  points  in  this  region  of  the  plot.  These 
inflections  generate  extrema  in  the  dispersion  of  group  velocities  and  will  be  associated  with  Airy 
phases  centered  on  frequencies  between  0.4  to  0.7  Hz.  An  examination  of  the  phase  velocities, 
group  velocities,  and  ellipticity  of  the  Rayleigh  inodes  between  these  inflection  points  suggests  that 
the  inflections  are  most  likely  to  be  products  of  P  to  SV  coupling  and  conversion  within  the  steep 
suificial  gradient  layer.  Since  the  Langer  approximation  does  not  have  these  inflections  in  this 
region  of  the  dispersion  plot,  it  will  be  unable  to  account  for  the  frequency  dependent  coupling 
and  conversion  of  P  to  SV  and  SV  to  P  waves  in  the  gradient  layer.  This  is  primarily  a  problem, 
however,  at  large  values  of  X/L.  Comparisons  of  Rayleigh  and  Love  mode  dispersion  for  models  2 
and  3  shown  in  Figure  2  found  no  such  discrepancies  when  X/L  was  less  than  0.1. 

Scalar  Amplitude  Function  and  Eigenfunctions 

The  scalar  amplitude  functions,  RA,  lA ,  which  are  generated  by  a  residue  calculation  at  a  pole 
in  the  complex  wavenumber  plane,  are  related  to  a  depth  integral  of  surface  normalized  energy 
(llarkrider  and  Anderson,  1966),  i.e.,  RA  =  f™  p[RE\(z)  +  REl(z)\/ rEx(G)  dz.  Figure  4 
compares  the  scalar  amplitude  function,  RA,  for  the  fundamental  Rayleigh  mode  in  models  1,  2, 
and  3  calculated  using  the  Langer  approximation  in  the  surface  gradient  layer  ( dashed  line )  versus 
those  calculated  by  parameterizing  the  gradient  layer  by  thin  layers  ( solid  line).  Once  again,  it  is 
clear  that  errors  in  the  use  of  the  Langer  approximation  become  smaller  as  X/L  decreases.  Errors 
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remain  sufficiently  small  for  practical  calculations  for  XfL  less  than  0.4.  Errors  for  model  3,  having 
the  lowest  magnitude  of  gradients  in  the  surface  layer,  remain  small  at  much  higher  values  of 
XfL  than  those  errors  for  models  1  and  2.  This  relative  behavior  can  be  explained  by  the  energy 
distribution  of  the  fundamental  mode  in  the  different  models.  The  energy  of  the  fundamental 
mode  is  confined  to  the  gradient  layer  in  model  1  throughout  the  plotted  wavelength  band,  whereas 
progressively  more  energy  is  distributed  in  the  zero  gradient  halfspace  in  models  2  and  3.  Hence, 
for  the  calculation  of  the  scalar  amplitude  factors,  the  breakdown  in  the  Langer  approximation 
depends  not  only  the  size  of  scale  length  of  a  gradient  layer,  but  also  on  the  distribution  of  energy 
with  depth  of  specific  modes  in  layers  having  strong  gradients. 

The  error  in  the  eigenfunctions  of  individual  modes  (e.g.,  Figure  5)  becomes  acceptably  small 
(less  than  5%)  when  X/L  less  than  0.4.  The  results  of  these  tests  can  be  summarized  in  terms 
of  frequency  by  a  rigorous  upper  bound  of  5%  error  in  computed  pole  positions,  scalar  amplitude 
functions,  and  eigenfunctions  for  frequencies  /  >  5  |W|. 

Synthetic  Seismograms 

Seismograms  are  synthesized  by  summing  modes  over  a  broad  band  of  frequencies.  An  upper  bound 
of  error  in  a  synthetic  seismogram  can  thus  be  based  on  the  tests  of  errors  of  pole  position,  residue 
functions,  and  eigenfunctions  using  either  the  lowest  frequency  (longest  wavelength)  in  a  broad 
band  synthetic  or  the  dominant  frequency  (wavelength)  of  a  narrow  band  synthetic. 

Figures  6  and  7  compare  reference  synthetics  and  Langer  approximated  synthetics  for  models  1 
and  3  shown  in  Figure  2.  Following  the  examples  shown  in  Spudich  and  Ascher  (1983),  the  seismo¬ 
grams  were  synthesized  to  represent  the  far-field  particle  velocity  of  a  step  function  in  displacement 
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convolved  with  a  low  pass  filter  that  rolls  off  with  a  cosine  taper  at  0.5  Hz.  Thus,  the  pass  band  of 
the  synthetic  seismograms  has  a  peak  centered  at  about  0.5  Hz. 

Although  kinematic  errors  in  the  mode  dispersion  calculations  arc  generally  small  throughout 
most  of  the  frequency  band,  the  dynamic  errors  in  mode  amplitudes  in  model  1  are  sufficient  to 
produce  poor  matches  in  the  group  velocity  band  corresponding  to  the  fundamental  mode  and  the 
first  few  higher  modes.  These  effects  can  be  seen  in  Figure  6,  in  which  the  early  portion  of  the 
seismograms  computed  by  the  two  methods  are  more  closely  in  phase  but  become  progressively 
out  of  phase  in  the  time  window  corresponding  to  the  arrival  of  the  fundamental  mode  and  first 
few  higher  modes.  The  agreement  between  tne  two  methods  is  much  better  for  the  transverse 
component  than  for  the  radial  or  vertical  components  of  motion,  undoubtedly  because  P  to  S 
coupling  in  the  gradient  layer  can  be  ignored  in  the  Love  mode  synthesis. 

The  results  (Figure  7)  for  model  3,  which  has  the  weakest  gradient,  show  the  match  between  ref¬ 
erence  and  Langer  approximated  synthetics  becoming  nearly  perfect.  Since  seismograms  computed 
by  the  two  methods  overlay  one  another  to  within  the  thickness  of  plotted  lines,  any  differences 
can  be  highlighted  only  by  plotting  the  difference  between  the  two  synthetics.  The  difference  seis¬ 
mograms  tend  to  be  largest  whenever  there  are  time  shifts  in  peak  oscillations.  Very  small  time 
shifts  (smaller  than  the  thickness  of  plotted  lines)  can  produce  visible  bumps  in  the  difference  seis¬ 
mograms.  At  the  dominant  frequency  of  the  synthetic  seismograms,  X/L  equals  0.4.  Although  the 
differences  between  the  synthetics  are  quite  acceptable  in  an  overlay  plot,  a  reduction  of  X/L  by  one 
half  or  more  would  be  required  to  nroduce  acceptably  flat  difference  seismograms.  A  conservative 
conclusion  is  that  errors  in  the  use  of  the  Langer  approximation  become  less  than  5%  when  the 
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ratio  \/L  is  less  than  or  equal  to  0.2. 


INTRINSIC  ATTENUATION 


To  be  practically  useful,  any  method  of  synthesizing  complete  seismograms  at  local  and  regional 
distances  must  be  capable  of  including  intrinsic  attenuation.  The  incorporation  of  the  attenuation 


in  the  Langer  approximation  simply  consists  of  the  analytic  continuation  of  all  formulae  to  complex 
velocities  (Cormier  and  Richards,  1976,  1989).  Care  must  be  exercised  in  the  definition  of  branch 
cuts  of  square  roots  and  fractional  powers  appearing  in  both  the  analytic  expressions  and  function 
subroutines  used  in  evaluating  the  Langer  approximation  (see  Appendix  I),  but  this  is  not  an 
insurmountable  problem.  The  Langer  subroutine  modified  for  use  with  locked  mode  calculations 
has  been  tested  in  problems  involving  integration  in  the  complex  ray  parameter  plane  combined  with 
complex,  frequency  dependent  velocity.  It  returns  generalized  vertical  wavefunctions  and  slownesses 
that  are  continuous  in  the  complex  ray  parameter  plane  except  for  poles  and  branch  cuts,  which 
emanate  from  complex  ray  parameters  corresponding  to  grazing  incidence  on  boundaries  in  an 
anelastic  model.  Test  calculations  have  demonstrated  that  the  position  of  these  singularities  does 
not  impede  a  successful  search  for  the  complex  zeros  of  the  dispersion  functions  of  locked  modes  in 
an  anelastic  model. 
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A„  = 

\tu  +  OJ2J 

Ar  .  ,n(‘j*±!l) 

\tur  +  u2y 


(5) 


(3r  is  the  real  velocity  at  a  reference  frequency  ur.  Equation  5  is  uniformly  valid  both  in  the 
center  of  the  absorption  band  as  well  as  across  and  through  the  high  and  low  frequency  corners  of 
the  absorption  band.  Complex  P  velocity  a  is  calculated  by  the  same  formula,  with  an  option  to 
constrain  attenuation  to  be  pure  shear  or  to  specify  a  different  peak  attenuation  parameter  Q~ 1  for 
P  waves.  Ideally  the  reference  frequency  uT  should  be  chosen  to  be  in  the  middle  of  the  frequency 


band  of  the  seismic  data  used  in  determining  a  trial  model  for  a  given  region.  Complex  velocities  are 
calculated  at  each  layer  boundary  by  equation  5  above,  and  linear  gradients  of  complex  velocity 
are  assumed  in  each  layer.  The  delay  time  function  r  needed  by  the  Langer  approximation  is 
calculated  as  described  in  Appendix  II,  but  it  now  must  be  recalculated  at  each  frequency.  It  is 
possible  to  specify  different  peak  Qp  values  as  well  as  different  upper  and  lower  limits,  u\  and  o>2, 
of  the  relaxation  band  at  the  top  and  bottom  of  each  inhomogeneous  layer. 

A  test  anelastic  model  is  shown  in  Figure  8.  The  attenuation  model  is  an  absorption  band 
model  in  pure  sh^ar  attenuation  having  gradients  in  peak  attenuation  Q^1,  and  low  and  high 
frequency  corners,  u>i,  u)2,  of  the  relaxation  band.  A  minimum  value  of  Qp  =  20  is  assumed  at  the 
surface.  The  velocities  and  Q  values  are  similar  to  values  measured  from  regional  seismograms  in 
New  England  (Kafka  and  Reiter,  1987).  Locked  mode  seismograms  were  synthesized  in  this  model 
using  two  different  approaches.  In  the  first  approach,  only  the  real  part  of  the  complex  velocities 
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is  used  in  calculation  of  mode  amplitudes  and  eigenfunctions,  and  a  complex  phase  velocity  was 
substituted  in  the  cylindrical  harmonics  describing  the  horizontal  propagation  of  each  mode.  This 
complex  phase  velocity  is  taken  from  the  complex  pole  k  estimated  by  first  order  perturbation 
theory.  This  is  the  standard  approach  for  handling  attenuation  in  surface  wave  and  locked  mode 
calculations  (Harvey,  1985;  Panza  and  Sudhadolc,  1987),  and  is  assumed  to  be  accurate  if  the  Q 
factor  is  sufficiently  high.  Day  et  al.  (1989)  have  shown  this  approach  to  be  inaccurate  for  some 
regional  seismic  phases  even  at  Q  values  on  the  order  of  several  hundred.  For  this  reason,  an 
exact  approach  was  developed,  in  which  a  search  was  made  for  the  complex  roots  of  the  dispersion 
function  and  all  formulae,  including  amplitude  factors  and  eigenfunctions,  were  evaluated  at  these 
complex  roois.  The  complex  pole  searching  algorithm  was  based  on  one  suggested  by  Schwab  and 
Knopoff  (1971),  with  modifications  near  osculating  points  of  the  dispersion  curves.  Near  these 
points,  the  complex  roots  are  found  by  the  same  algorithm  for  a  series  of  increasing  Q'1  values, 
approaching  the  true  Q~ 1  model.  Checks  are  made  for  duplication  or  omission  of  poles  at  the  end 
of  this  procedure  for  each  frequency. 

Figure  9  compares  the  results  of  these  two  methods  for  incorporating  attenuation  of  the  funda¬ 
mental  mode  Rayleigh  wave.  The  seismograms  computed  by  the  different  methods  nearly  overlay 
one  another  at  all  distances.  The  exact  method  reduces  some  high  frequency  numerical  noise, 
which  is  barely  visible  at  the  scale  of  Figure  9.  The  differences  in  the  complex  phase  velocities 
computed  by  the  two  methods  are  on  the  order  of  0.001  km/sec  in  the  real  part  of  the  complex 
phase  velocity  and  vary  from  1  x  lO"10  to  1  x  10-4  km/sec  in  the  complex  part  of  phase  velocity 
as  frequency  increases  up  to  2  IIz.  The  differences  between  the  depth  behavior  of  the  real  part  of 
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the  complex  eigenfunctions  are  insignificant  between  the  two  methods.  From  these  results  it  can 
be  concluded  that  the  perturbation  approach  to  attenuation  remains  very  accurate  in  the  synthesis 
of  the  fundamer  ta!  mode  for  Q  values  as  low  as  20.  For  the  synthesis  of  higher  modes,  particu¬ 
larly  those  contributing  to  refracted  P  and  S  and  interference  head  waves,  more  detailed  tests  have 
shown  that  the  perturbation  approach  introduces  significant  error  as  Q  values  decrease  below  100. 
Experiments  in  progress  by  an  author  of  this  paper  (Harvey)  document  a  progressive  breakdown  in 
the  perturbation  approach  as  phase  velocity  increases,  corresponding  to  more  vertically  propagat¬ 
ing  waves  in  attenuating  layers.  The  safest  approach  to  handle  attenuation  accurately  in  a  locked 
mode  method  is  to  develop  all  calculations  around  the  exact  complex  modes.  Our  motivation  in 
discussing  this  example  is  simply  to  show  that  an  exact  complex  mode  search  can  succeed  with  the 
use  of  the  Langer  approximation  together  with  complex  dispersive  velocities.  Use  of  the  Langer 
approximation  easily  permits  the  implementation  of  the  assumption  that  gradients  in  the  real  part 
of  the  elastic  moduli  are  also  associated  with  gradients  in  the  imaginary  part  of  elastic  moduli. 

It  is  appropriate  to  conclude  this  section  with  some  cautionary  words.  Often  a  very  low  Q 
layer  is  required  in  a  surface  layer  in  order  to  produce  realistic  simulations  of  seismograms  observed 
at  local  and  regional  distances  (e.g.,  Panza  and  Sudhadolc,  1987).  If  the  apparent  attenuation 
of  such  a  layer  is  truly  due  to  viscoelasticity,  its  effects  can  be  accurately  calculated  by  complex 
locked  modes.  It  is  more  likely,  however,  that  such  apparent  low  Q's  are  due  to  a  combination  of 
scattering  by  topography  of  layer  boundaries  and  volumetric  heterogeneities  and  frictional  sliding 
of  grains  and  open  cracks.  Neit.’.-ir  of  these  effects  can  be  simulated  by  a  combination  of  vertically 
varying  layers  and  linear  viscoelastic  relaxations. 
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EFFECTS  OF  GRADIENTS  ON  REGIONAL  PHASES 


To  test  the  effects  of  crustal  and  upper  mantle  gradients  on  regional  seismic  phases,  locked  mode 
synthetics  were  computed  in  two  simple  models  MH  and  MG  (Figure  10).  Model  MH  consists 
of  a  two-layered  crust  overlying  a  homogeneous  mantle.  MH  has  also  been  used  for  testing  and 
benchmark  timing  of  many  different  techniques  of  computing  complete  seismograms  at  regional 
distances  (Richards  and  Mithal,  personal  communications).  Model  MG  consists  of  a  single  crustal 
layer  having  a  positive  gradient  with  depth,  overlying  a  mantle  having  a  positive  gradient  with 
depth.  The  mantle  gradient  is  consistent  with  the  increase  in  seismic  velocities  typical  of  reference 
earth  models  between  the  Moho  and  400  km  depth.  The  depth  averaged  crustal  velocities  of  MH 
and  MG  are  identical.  Both  models  have  an  attenuation  structure,  with  Q’s  in  a  high  enough  range 
that  simple  perturbation  theory  can  be  accurately  used  to  calculate  the  effects  of  attenuation  in 
the  locked  mode  method.  Seismograms  were  synthesized  in  a  frequency  band  up  to  2  Hz  for  the 
source  and  receiver  geometries  used  by  W-Y.  Kim  (1987),  who  synthesized  seismograms  in  model 
MH  using  wavenumber  integration. 

The  synthetic  seismograms  (Figure  11)  for  the  first  10  higher  Rayleigh  modes  (fundamental 
mode  excluded  from  sum)  are  very  similar  in  peak  amplitude  and  length  of  coda  in  both  models 
MH  and  MG.  The  group  velocity  window  of  the  energy  centroid  corresponds  to  that  expected  for  the 
Lg  phase.  The  strong  similarity  of  the  synthetic  seismograms  suggests  that  Lg  is  not  very  sensitive 
to  the  details  of  the  crustal  model,  its  coda  primarily  being  controlled  by  the  total  thickness  of  the 
crust  and  its  average  shear  velocity.  It  is  probably  possible  to  simulate  realistic  Lg  phases  using  a 
very  few  number  of  crustal  layers.  Introduction  of  crustal  layers  in  a  modeling  experiment  may  not 
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be  necessary  unless  there  is  compelling  evidence  for  crustal  discontinuities  observed  in  the  earlier 
time  window  in  the  form  of  refracted  body  waves  and  interference  head  waves. 

To  illustrate  the  effects  of  gradients  in  the  mantle  and  to  provide  a  further  check  on  the  accu¬ 
racy  of  the  Langer  approximation  in  more  complicated  models  at  longer  range,  seismograms  were 
synthesized  in  models  MU  and  MG  using  both  the  Langer-locked  mode  method  and  a  wavenumber 
integration  method  (Mandal  and  Mitchell,  1986)  in  which  gradient  layers  were  parameterized  by 
thin  homogeneous  layers.  For  the  wavenumber  integration  method,  about  100  thin  layers  were  used 
to  compare  with  the  results  of  the  Langer-lc  ~ked  mode  method  (1.5,  3,  and  7  km  thick  layers  were 
used  for  the  upper  38  km,  38-128  km,  and  deeper  than  128  km,  respectively).  In  this  frequency 
band,  ringing  effects  were  observed  in  the  wa  'enumber  integration  method  when  the  thin  layers 
simulating  the  gradient  zones  exceeded  10  km  in  thickness. 

The  close  agreement  (Figures  12  and  13)  in  the  synthetics  calculated  using  quite  different 

approaches  and  model  parameterizations  demonstrate  that  the  Langer  approximation  is  sufficiently 

accurate  for  the  gradients  of  model  MG.  The  maximum  magnitude  of  gradient  in  model  MG  occurs 

in  the  P  velocity  of  the  mantle,  which  is  0.0038  sec"1.  At  this  gradient,  errors  in  the  Langer 

approximation  are  bounded  by  5%  for  frequencies  /  >  0.016  Hz  or  wavelengths  A  <  480  km.  The 

lowest  observable  frequency  component  in  the  pass  band  shown  in  Figures  11-13  is  about  0.2  Hz. 

Any  differences  between  the  synthetics  for  model  MH  and  MG  are  thus  truly  due  to  differences 

in  wave  propagation  in  the  two  models  and  not  to  the  method  of  synthesis.  We  found  the  peak 
amplitudes  obtained  with  the  plane  layered  method  to  be  extremely  sensitive  to  the  thickness  of 
the  thin  layers  chosen  to  simulate  the  gradient  layers.  Thus  we  attribute  small  differences  in  the 
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peak  amplitude  scaling  of  the  two  methods  to  interference  effects  of  thin  layering.  In  a  comparison 
of  the  results  for  models  Mil  (Figuie  12)  and  MG  (Figure  13),  it  is  seen  that  the  seismograms  are 
very  similar  at  closer  ranges  but  at  300  km  some  differences  begin  to  be  notable.  Pn  and  Sn  are 
very  weak  in  the  Mil  simulation,  but  are  very  strong  in  the  MG  synthetic.  Pn,  Sn,  and  crustal 
reverberations  converted  to  Pn  and  Sn  are  so  strong  in  the  MG  synthetic  that  they  dominate  Lg  in 
amplitude,  and  the  seismogram  becomes  a  series  of  spikes  consisting  of  P  and  S  reverberations  in 
the  crust  converted  to  Pn  and  Sn  interference  head  waves  in  the  mantle.  The  comparison  confirms 
what  is  known  about  Pn  and  Sn  as  interference  head  waves  in  models  having  positive  gradients 
below  the  Moho.  A  positive  gradient  acts  to  enhance  the  amplitude  of  the  interference  head  wave 
far  above  what  would  be  predicted  for  a  classical  head  wave  in  a  homogeneous  layer  (Hill,  1971; 
Cerveny  and  Ravindra,  197T,  Menke  and  Richards,  1980).  Although  the  steep  mantle  gradients  of 
model  MG  may  not  be  appropriate  for  many  regions,  they  are  within  the  probable  range  of  lateral 
variations  in  upper  mantle  structure.  As  noted  by  Hill  (1971),  a  tradeoff  exists  between  Q  and 
upper  mantle  velocity  gradients  in  the  modeling  the  amplitude  of  Pn  and  Sn.  For  example,  the 
relative  excitation  of  Pn  and  Sn  versus  Lg  in  model  MG  can  be  made  to  appear  more  similar  to 
that  in  model  MH  by  decreasing  Q  in  the  mantle. 

CONCLUSIONS 

In  this  paper,  we  have  demonstrated  that  even  small  gradients  of  VV  =  0.004  sec-1  can  substan¬ 
tially  affect  the  distance  decay  of  interference  head  waves  such  as  Pn  and  Sn.  Lg,  on  the  other 
hand,  is  only  verly  weakly  sensitive  to  details  of  crustal  layering  or  gradients.  The  peak  amplitude 
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positions  and  eigenfunctions  compared  to  those  synthesized  using  the  Langer  approximation  in  the 
gradient  layer  parameterized  by  analytic  velocity  functions.  A  calculation  in  a  thick  homogeneous 
layer,  however,  would  always  be  more  efficient  than  a  calculation  using  the  Langer  approximation 
in  an  inhomogeneous  layer  of  the  same  thickness.  A  model  parameterization  that  may  be  the  best 
compromise  between  computational  efficiency  and  realism  in  the  behavior  of  regional  phases  would 
be  one  having  a  crust  composed  of  homogeneous  layers  overlying  a  mantle  composed  of  gradient 
layers.  Seismograms  synthesized  in  such  a  model  could  accurately  predict  the  Lg  phase  as  well  as 
the  Pn  and  Sn  phases.  (This  study  did  not  investigate  the  importance  of  crustal  gradients  for  the 
Pg  phase.) 

The  Langer-locked  mode  approach  to  synthesizing  complete  seismograms  may  also  offer  some 
advantages  in  waveform  inversion  for  earth  structure.  By  reducing  the  number  of  parameters  needed 
to  describe  a  model,  the  inverse  problem  for  structure  would  be  simplified  and  fewer  experiments 
would  be  needed  to  determine  the  maximum  number  of  resolvable  layers.  A  layer  need  only  be 
introduced  when  the  data  firmly  suggest  the  existence  of  first  order  discontinuities.  Furthermore, 
specification  of  the  earth  model  by  a  small  number  of  gradient  layers,  bounded  by  well  known  first 
or  second  order  discontinuities,  facilitates  a  comparison  of  inverted  parameters  with  the  velocity, 
attenuation,  and  density  behavior  predicted  by  theoretical  or  empirical  formulae  given  as  functions 
of  pressure  and  temperature. 
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APPENDIX  I  -  THE  L ANGER  APPROXIM  A  TION 


Vertical  Wavefunctions 

Tlie  notation  for  the  Langer  approximation  (Langer,  1932;  1949)  differs  among  different  authors 
who  have  applied  it  to  seismic  wave  propagation.  (Richards,  1976;  Woodhouse,  1978;  Chapman, 
1974;  Doornbos,  1981),  involving  either  llankel  functions  of  order  1/3  or  Airy  functions  of  different 
types  or  arguments  to  give  exponentially  decaying  and  growing  type  solutions  below  a  turning 
point.  The  notation  adopted  here  is  basically  that  given  in  Aki  and  Richards  (1980). 

The  Langer  approximation  is  a  uniformly  asymptotic  approximation  to  the  vertically  separated 
part  of  the  solution  to  the  elastodynamic  wave  equation  in  a  region  in  which  elastic  moduli  and 
density  vary  continuously  with  depth.  The  zeroth  order  term  in  frequency  in  the  asymptotic  solution 
is  given  as 


jr^(r)  =  Vine  3 


. .,  2 is 
Ai{-  e  3 


GO 


ff(2)(r)  =  y/2n  e’3 


*  1  or 

X'\ 

4/4' 


Ai(-  e  3  (a) 

*(3)W  =  ^  (££)  M-(.) 


a(1)(r)  =  v/Tjr  eT1 
CT^(r)  =  \/2  7r  cT 

CT^(r)  =  spun 


($) 


(AI.l) 
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where 


At  is  an  Airy  function  and 


<«  =  (3/2  ura)2'3 


^  =  (3/2  urtf'3 


Xa  =  sjl/a2  -  Jp/r2 


Xp  =  yJl/P2  -  f/r2 


a  and  P  are  the  P  and  S  velocity  respectively  at  radius  r,  p  is  the  ray  parameter  in  a  spherical 
Earth,  and  rp  is  the  turning  point  radius,  i.e.,  that  radius  at  which  Xa  or  A^  vanishes.  In  each 
inhomogeneous  layer,  the  velocity  functions  a(r)  and  (3(r)  are  assumed  to  be  analytic  and  to  possess 
only  one  turning  point  rp  in  the  domain  of  complex  p  used  in  synthesizing  a  seismogram. 

The  7r  wavefunctions  are  those  for  P  waves;  the  cr  wavefunctions  are  those  for  S  waves.  Several 
possible  pairs  of  independent  solutions  may  be  chosen  to  define  fundamental  matrices,  which  can 
be  used  to  solve  problems  in  wave  propaga*ion  in  media  consisting  of  a  sequence  of  vertically 
inhomogeneous  layers.  The  pairs  (tt^,  jt^)  and  (crM,  correspond  to  up-  (1)  and  down¬ 
going  (2)  waves.  The  pairs  (?r^,  x^)  and  (cr^\  a^)  correspond  to  up-going  (1)  and  standing 
or  evanescent  waves  (3).  When  the  turning  point  radius  rp  is  greater  than  r  and  ^r(^)  is  positive 
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with  decreasing  r,  the  r,  the  wavefunctions  irt3)  and  arc  exponentially  decaying  functions  with 
decreasing  r. 


Vertical  Slownesses 


Implementation  of  the  Langer  approximation  in  problems  in  which  clastic  boundary  conditions 
must  to  be  satisfied  at  model  discontinuities  is  simplified  by  the  introduction  of  generalized  cosines 
(Richards,  1976;  Aki  and  Richards,  1980)  or  generalized  vertical  slowness  functions,  which  are 
defined  as  follows 


\ 

i 


dir^1) 

dr 

dirt2) 

dr 

_  dirt 3) 
dr 


/(iuiirt1)) 
/( iwirt h) 

/(iuirt3)) 


(A1.2) 


do^  ... 


v  =  - 
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dr 
dot2) 
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dr 


/(iuot1)) 

/(iuot3)) 


The  normalization  of  the  vertical  wavefunctions  differs  slightly  from  that  given  in  Aki  and 
Richards  (1980)  and  has  been  chosen  such  that  the  following  relations  are  satisfied 


4*  ^7T^7T^  =  1 

£ir  f1)^3)  +  i-irt3)  x(i)  —  i 
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(AI.3) 

*r<»>*<a>  +  r^M1*  =  1 
^(l)ff{3)  +  ^(3)ff(l)  =  j 

These  relations  can  be  demonstrated  by  substituting  the  Langer  approximation  for  the  vertical 
slownesses  and  using  the  Wronskian  relations  between  the  Airy  functions  having  different  argu¬ 
ments.  Equations  AI.3  are  satisfied  exactly  when  only  the  zero  order  terms  in  frequency  are  kept 
in  the  definitions  of  the  vertical  slownesses. 

Fundamental  Matrices 

Boundary  conditions  in  a  medium  consisting  of  n  inhomogeneous  layers  can  be  handled  in  the  same 
manner  as  a  medium  consisting  of  homogeneous  layers,  but  with  the  Langer  approximation  to  the 
vertical  wavefunctions  and  vertical  slownesses  substituting  for  exponential  functions  and  cosines. 

P-SV 

As  a  function  of  radius  r,  the  fundamental  matrix  for  P-SV  propagation  and  Rayleigh  modes  is 
taken  to  be  that  given  in  Cormier  (1980): 

(fjrM  p/r  oW  p/roW 

—ip/rirW  —ip/rirW  ifjaM  —iijcrW 
-t'Ajrl1)  -iAirW  iBrp(l)  -iBr\o{2) 

Bi*M  -B$1 r(2>  AaW  Aa& 

(AI.4) 
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iAx^  B£pffl  —  —  i  p/r 

-iAn^  ip/rv M 

-iBijcrW  AoW  —  p/r  ifja^ 

-iBr/oW  —AoW  p/r  <r(l)  ifp^ 

The  fundamental  matrix  may  alternatively  be  defined  using  the  wavefunction  pairs  (7^),  v^) 
and  (Cormier,  1980).  This  fundamental  matrix  has  exactly  the  same  form  as  AI.4, 

but  with  (3)  replacing  the  (2)  superscripted  wavefunctions  and  the '  accent  replacing  the '  in  the 
vertical  slownesses.  In  all  calculations,  the  (3)  superscripted  wavefunctions  are  substituted  for  the 
(2)  superscripted  (down-going)  wavefunctions  in  the  p  domains  in  which  exponentially  decaying 
and  growing  vertical  wavefunctions  exist.  With  a  few  simple  modifications  described  by  Doornbos 
(1981),  the  fundamental  matrix  defined  in  AI.4  can  be  applied  to  layers  having  a  negative  as  well 
as  a  positive  gradient  with  depth. 


Fundamental  Matrix  for  SH  Propagation 


The  SII  fundamental  matrix  and  its  inverse  are 


F(r)  =  V 


p-W  <r(2) 
-*y/2  fjai 2) 


F(r)-1  =  V 


-i  /i1/2  r/< r(2)  -/i-1/2  or( 2) 

y/V1) 


(AI.5) 
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Model  Parameterization 


Since  the  Langer  approximation  allows  layers  to  be  vertically  inhomogeneous,  the  effects  of  Earth 
curvature  are  built  into  the  model  parameterization.  All  formulae  are  evaluated  using  velocities 
and  densities  given  as  functions  of  radius,  r,  from  the  Earth’s  center.  In  each  inhomogeneous  layer, 
the  velocities  are  specified  by  analytic  functions,  which  have  only  one  turning  point  solution  in  the 
p  domain  of  interest.  Layer  boundaries  are  introduced  and  boundary  conditions  are  evaluated  at 
discontinuities  in  velocity  derivatives  as  well  as  first  order  discontinuities. 

To  provide  analytic  forms  for  the  delay  time  functions  ra  and  rp,  each  inhomogeneous  layer  is 
parameterized  by  making  the  flattened  velocity  be  a  linear  function  in  the  flattened  depth  coordi¬ 
nate,  z.  The  usual  (Muller,  1971)  mapping  between  the  flattened  velocity  function  vj(z)  and  the 
true  velocity  function  u(r)  is  assumed  : 


«(r)  =  rvj{z)/Re 

where 

z~  =  Re\og(rc/Re) 
where  Re  is  the  radius  of  the  Earth. 

The  flattened  velocity  function  vj  is  assumed  to  be  a  linear  function  in  flattened  depth,  com¬ 
puted  from  the  values  of  vj  at  flattened  depths  z^  and  z+_l  corresponding  to  radii  r~  and  lt 
bounding  the  top  and  bottom,  respectively,  of  vertically  inhomogeneous  layer  n.  The  analytic  form 
of  the  delay  time  function  r(r)  becomes 
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two  depths  such  that  analytic  integration  of  the  r  integral  becomes  possible.  Alternative  model 
parameterizations,  which  give  an  analytic  form  of  r,  are  discussed  by  Cormier  (1980),  Cerveny  and 
Jansky  (1983),  and  Cormier  and  Richards  (1989). 
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APPENDIX  H  -  MODE  AMPLITUDES  AND  EIGENFUNCTIONS 


Rayleigh  Modes 

The  summation  of  locked  Rayleigh  modes  requires  the  calculation  of  an  antisymmetric  Y  matrix 
having  five  independent  elements. 

The  Y  Matrix 

At  the  radius  rc  at  the  top  of  the  capping  layer,  starting  values  of  the  Y  matrix  are  taken  to  be 

Yu  =  —  A2  —  B2  Xac  Xpc 

V13  ~  -Ac  p/re  —  Bc  Xac  Xpc 

Yu  =  -i  Pc  Xpc  (AII.l) 

1 23  =  *  Pc  Xac 

V34  =  -  Xpc  Xac  -  p2/rl 

where  t  =  V^T,  and 


V.  =  *  sjp2/rl  -  1/P2 

K  =  *  \Jp2/rc  -  i/qJ 

Ac  =  2p2/r]  nc  -  Pc 


Bc  =  2  p2/r2c  pc 
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and  /zc,  pc,  <*c,  Pc  are  the  shear  modulus,  density,  P  velocity,  and  S  velocity,  respectively,  of  the 
high  velocity  capping  layer.  With  the  Yij  defined  here,  the  scalar  amplitude  function  rA  is  related 
to  that  defined  in  Harvey  (1981),  rAh,  by  rAh  =  -  n\p/Re. 

At  any  radius  r,  Y  can  be  computed  from  the  product 

Y(r)  =  KT(r,r+)Y(r;)K(r,r+)  (AII.2) 

where  K  is  a  P-SV  propagator  matrix  equal  to  a  product  of  intralayer  propagator  matrices  for  each 
layer,  m,  m  +  1,  etc. 

K  =  Km(r,r+)  Km+1(r~,r++1) - Kn(r;_„r+)  (AII.3) 

Layers  are  separated  by  boundaries  at  which  velocities  and/or  densities  have  either  first  or 
second  order  discontinuities.  Within  each  layer,  the  velocity  functions  are  continuous,  analytic 
functions.  Each  interlayer  propagator  matrix,  Km  is  constructed  from  the  zeroth  order  term  in 
frequency  of  the  uniform  asymptotic  approximation  to  the  fundamental  matrix  F  of  the  inhomo¬ 
geneous  layer.  Since  the  uniform  asymptotic  approximation  of  Langer  is  assumed,  the  velocity 
functions  within  each  layer  must  have  no  more  than  one  turning  point  for  each  ray  parameter,  p. 
With  this  restriction,  computations  can  still  be  conducted  in  a  complicated  model  having  one  or 
more  low  velocity  zones,  as  long  as  this  model  is  built  from  ’’layers”  in  which  the  analytic  functions 
for  P  and  S  velocity  have  only  a  single  turning  point  for  each  p. 

The  intralayer  propagator  is  defined  by 

K.K.4J  =  ^(r-)  F-1(r++1)  (AIM) 

Substituting  in  equation  AII.2  the  forms  for  the  fundamental  matrix  and  its  inverse  from 
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equation  AI.4,  and  simplifying  the  resulting  expression  gives  recursion  relations  as  follows  for  the 
upward  propagation  of  Y  matrix  elements: 


Yi 2(0  =  £  kWn  kGn  -  2  A(r+_,)  J?(r+_,)  0Wn  0Gn 

k= 1 

=  »*(r;)  WrtiWrt,)  -  PlTt-\  «(r-i+)]*W'.»G, 

Jfesl 


^n1)  =  "E  hMtDkWnkQn 

k=l 


(AII.5) 


^(rn1)  =  -E  kd^r-1)  kWn  kGn 

k=  1 

Mr;)  =  -  E  *d»(rn ')  *Wn  +  2  p/f+Lj  O^n  O^n 

*=1 

where  the  quantities  kdi(z),  kW„,  kGn  are  defined  as  follows: 


kdi(r)  =  -A\r)  -  B\r)  kXa{r)  kX0 


kd2(r)  =  A(r)  p/r  +  B{z)  kXa(r)  kX „(r) 


kd3(r)  =  ip(r )  kX0(r) 


(AII.6) 


kd4(r)  =  —ip(r)  kXa(r) 
kds{r)  =k  A/3 (r)  *Aa(r)  +  (p/r)3 


kWn=  -*d,(r+_,)MV-i)  +  2*d2(r+_1)M»V-i)  + 


+  #ti)J,J3(Ci)  +  *rf5(r+.J)l'i2(r;.1)  (AII.7) 
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for  k  ^  0  and 


o Wn  =  p/r+  y,2(r;_,)  +  [A(r+_,)  +  p/r+. ,  B(r+^)]Yia( r~_,) 

rn/^-l 


oGn  =  4 


yjMrn  )pirt-i) 


i Gn  =  7r(2)(r;)  a(2)(r~)  ir(1)(r+_1)  <r(1)(r+_i)  / 

yjMrZWt-i) 


2<j„  =  -7T<2)(r„  )  cr(1)(r„)  7r(1)(r+_i)  <r(,)(r+_i)  y2^^+  ^ 


3Gn  =  -7r(1)(r")  <r<2)(r;)  ^(r+.j)  <r^>(r+_,)  -====== 


*Gn  =  JT(1)(rn)  ^(1)(rn)  ^(^-l)  ^(^-i)  /  ^"T*;  \ 

V2/>(rn)p(r+.i) 

*A0  and  jtA^  denote  the  following  at  boundaries  r~  and  r+  : 
lM^)  =  f(r“) 
iAa(r+.i)  =  £'(r+_i) 


l¥rn‘)  =  ^n) 

iA/3(r„_j)  =  v(rt-i) 

2  =  £(0 

2Aa(r+_1)  =  e'(r^-i) 
^pirn)  =  ~V(rn) 
2A/j(r+_i)  =  “4(rti) 

3^0  (r-)  =  -£(r") 


(AII.8) 


(AII.9) 
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3Aa(r+_,)  =  -£(r+_i) 
3b(rn)  =  fan) 

=  ’K'n-l) 
4*o,(r~)  =  -^(r") 

<A/j(rn)  =  -Tj{r-) 
<Mrn-i)  =  -v{r+_ ,) 


Layer  Reduction 

The  first  term  ( k  =  1)  in  the  summation  in  equation  AII.5  is  of  the  same  form  as  the  starting 
values  Y  matrix  in  the  capping  layer  in  regions  of  slowness  in  which  the  vertical  wavefunctions 
behave  exponentially.  When  this  first  term  is  exponentially  larger  by  several  orders  of  magnitude 
than  the  ( k  =  2, 3, 4, 5)  terms,  then  the  Y  matrix  calculation  may  be  started  at  a  higher  layer, 
taking  this  higher  layer  as  the  capping  layer.  This  procedure  of  layer  reduction  is  analogous  to  that 
described  in  homogeneously  layered  models  (Panza  and  Sudhadolc,  1987). 

The  Capping  Layer  -  To  Cap  or  Not  to  Cap? 

A  capping  layer  may  be  avoided  and  the  starting  values  of  the  Y  mitrix  defined  by  evaluating  all 
velocities  and  densities  in  AII.2  at  the  top  of  the  deepest  layer  (r  =  r~ )  and  setting  Xac  and  \pc 
to  the  Langer  generalized  vertical  slownesses  £(r“)  and  f/(r~),  respectively  (Cormier,  1980).  The 
generalized  vertical  slownesses  |  and  fj  contain  the  phase  information  needed  to  represent  turning 
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rays  in  the  deepest  layer.  In  this  choice,  the  complete  seismogram  can  be  synthesized  either  by 
( 1 )  numerically  integrating  along  a  contour  close  to  the  real  k  axis  or  by  (2)  summing  residues  of 
locked  modes  at  high  wavenumbers  and  integrating  along  the  real  k  axis  at  low  wavenumbers.  The 
poles  at  high  wavenumbers  represent  guided  waves  trapped  in  the  upper  layers,  while  the  integral 
along  the  real  k  axis  at  low  wavenumbers  represents  body  waves  reverberating  at  high  incidence 
angles  in  the  upper  layers,  leaking  into  and  bottoming  in  the  deepest  layer. 

In  many  problems,  the  seismogram  may  be  represented  entirely  by  locked  modes  because  the 
structure  provides  a  natural  capping  layer  if  all  waves  of  interest  bottom  above  or  are  totally 
reflected  by  some  natural  layer  at  some  depth.  In  the  natural  capping  layer,  the  limiting  forms  of 
the  generalized  slownesses  £  and  fj  at  largo  values  of  \ljt\  in  the  upper  half  complex  p  plane  are 
equal  to  the  vertical  slownesses  in  a  capping  layer,  A0c  and  A pc,  defined  in  AII.2.  Unfortunately, 
this  natural  capping  layer  can  be  as  deep  as  the  inner  core  boundary  for  body  waves  interacting 
with  upper  mantle  structure.  To  minimize  the  number  of  numerical  operations  in  the  calculation 
of  dispersion  functions  and  eigenfunctions  and  to  retain  the  simplicity  of  a  modal  representation, 
it  is  more  convenient  to  specify  a  capping  layer  having  artificially  high  velocities. 

Eigenfunctions 

Although  propagation  of  the  Y  matrix  elements  has  been  shown  to  be  numerically  stable  at  arbi¬ 
trarily  high  frequency  (Abo-Zena  1979;  Harvey,  1981),  numerical  problems  in  the  calculation  of  the 
Rayleigh  eigenfunctions  occur  if  E  is  calculated  by  multiplying  propagator  matrices.  One  approach 
to  this  problem  is  to  divide  a  layer  into  thin,  pseudo  layers,  and  rescale  the  propagator  matrix  after 
propagation  through  each  thin  layer.  Better  techniques,  however,  can  be  formulated,  which  do  not 
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require  the  introduction  of  additional  pseudo  layers. 

One  technique,  described  by  Harvey  (1985),  expresses  the  eigenfunctions  in  terms  of  Y  matrix 
elements  by  propagating  the  wavefield  upward  from  the  cap  layer  and  downward  from  the  free 
surface.  Thus,  since  the  calculation  of  Y  elements  is  numerically  stable,  so  is  the  calculation  of  the 
E  eigenfunctions.  In  this  technique,  eigenvalues  can  be  normalized  at  the  source  depth,  offering 
numerical  advantages  in  the  calculation  of  channel  waves  having  vanishingly  small  energy  outside 
of  a  waveguide. 

The  technique  used  here  also  does  not  require  pseudo  layering,  but  retains  the  standard  nor¬ 
malization  of  the  Ei  function  to  1  at  the  free  surface.  The  first  step  in  this  technique  is  to  recognize 
that  the  stress  eigenfunctions  E3  and  E4  can  be  calculated  from  the  displacement  eigenfunctions 
Ei  and  E2  by 


E3  =  -  Yu/ Ei  -  Yu/ Ym  E2 


Ea  =  =  Yi3/Ym  Ei  +  Yn/Yu  E2 

Using  these  relations,  the  four  equations  that  propagat?  the  E  vector, 


E(r)  =  K(r,r„)  E(r„) 


can  be  rewritten  as  two  equations  that  propagate  Ei  and  E2> 


E\(r) 

E2(r) 


L(z,*„) 


E\(rn) 


H  2) 


(AII.10) 


(AII.11) 


(All. 12) 
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and  the  two  equations  given  in  AII.10  between  the  displacement  eigenfunctions  and  stress  eigen¬ 
functions.  Using  the  Y{}  defined  in  this  paper,  the  eigenfunctions,  are  related  to  those  defined 


in  Harvey  (1981),  Ej1,  by  £,  =  E{r,  E2  =  -  3?, 


£3  =  -  R'Ez/p,  and  E4  =  ReE4  /p. 


A  new  2x2  propagator  matrix  L  is  defined  having  components 


L\\  =  1(\\  -  A"i3  yj 4/^34  +  A 14  ^13/^34 


1-12  =  K\2  -  K\z  i  24/^34  +  A"i4  Y2z/Yz4 

(AII.13) 

Z21  =  A'21  •  23  ^14/^34  +  A’24  Y\z/Yz\ 

£22  =  A'22  +  K23  Y24!  ^34  +  A'24  V23/I34 

To  ensure  numerical  precision  in  a  machine  calculation,  the  individual  propagator  elements  as 
well  as  the  recursion  formulae  in  AII.5  for  the  Y  matrix  elements  must  be  substituted  into  the 
definitions  of  the  the  L{j  clervients  in  AII.13,  a  fraction  formed  with  the  common  denominator  of 
I34,  and  the  numerator  of  the  fraction  simplified.  When  this  simplification  is  done,  it  is  seen  that 
all  numerator  terms  that  potentially  are  of  the  largest  exponential  order  cancel.  Although  many 
cancellations  occur,  the  resulting  expressions  for  the  L{j  elements  are  still  quite  lengthy  and  are 
not  given  here. 
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Love  Modes 


D\  and  Di 

In  this  case,  calculation  of  the  dispersion  function  D\  eigenfunction  vector  E  can  proceed  by  simple 
multiplication  of  propagator  matrices  without  I06S  of  numerical  precision.  The  vector  {D\,  D2) 
in  the  notation  of  Harvey  (1985)  is  equal  to  the  vector  Esh  in  the  notation  of  Cormier  (1980). 
In  the  capping  layer,  (D\t  D2)  is  simply  equal  to  the  first  row  of  the  inverse  fundamental  matrix 
for  SH  waves.  Any  constant  may  be  chosen  to  multiply  the  starting  value  of  (I?|,  D2),  since  this 
constant  will  cancel  in  the  definition  of  eigenfunctions  and  in  the  ratio  appearing  and 

the  expression  for  the  total  response.  Starting  values  of  Dj  and  D2  at  the  top  of  the  cap  layer  are 
thus  taken  as 


D\  =  -i  pc  0c  Vc 


(AII.14) 


d2  =  -i!  Pc 


D\  and  D2  are  propagated  upward  by  multiplication  of  SH  propagator  matrices.  Since  (£>1,  J92 ) 
are  related  to  the  inverse  fundamental  matrix,  one  must  right  multiply  the  starting  values  by  the 
SH  propagator  matrix. 


Di(r)  D\(rc) 

=  K(r,r+)  (AII.15) 

D2(r)  D2(rc ) 

With  the  Di  defined  here,  the  scalar  amplitude  function  iK  is  related  to  that  defined  in  Harvey 
(1981),  L\li,  by  lA"  =  LA p/Le. 


Eigenfunctions 


Love  wave  eigenfunctions  are  defined  by 


£i(r)  =  D2(r)/D2(Re ) 


(AII.16'1 


E2(r)  =  D\(r)l  D2(Re) 

Using  the  D,  above,  the  eigenfunctions,  are  relateu  to  those  defined  in  Harvey  (1981),  E-1, 
hy  E{  =  Ei  and  E2  -  ReE2  /p.  In  the  residue  calculation,  scale  factors  can  be  applied  in 
each  layer  and  discarded  during  upward  propagation.  This  is  because  all  st  ale  factors  cancel  when 
rati0  m$h  is  formed.  In  the  eigenfunction  calculation,  the  total  scale  factor  of  each  A  must 
be  saved  in  order  to  describe  properly  regions  of  exponential  decay  of  the  eigenfunction.  In  the 
cases  where  E\  and  E2  are  exponentially  small,  the  depth  of  the  capping  layer  can  be  raised  and 
calculations  started  at  a  shallower  depth. 


Branch  Cuts  and  Poles  in  the  Complex  k  Plane 

The  functions  that  define  the  generalized  vertical  wavefunctions  and  slownesses  contain  branch  cuts 
emanating  from  points  in  the  complex  p  plane  corresponding  to  ray  parameters  grazing  the  model 
discontinuities.  Extreme  care  must  be  exercised  both  in  the  definition  and  the  choice  of  branch 
cuts  appearing  in  all  funeU  ms  oi  variables  raised  to  fractional  powers.  A  subroutine  has  been 
designed  such  that  substitution  of  Langer  generalized  wavefunctions  and  slownesses  into  coefficients 
of  refl-  non,  transmission,  and  conversion  always  gives  expressions  that  are  analytic  everywhere 
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in  the  complex  p  or  k  plan  „■  except  for  discrete  poles  and  zeros.  This  Langer  subroutine  has  been 
tested  in  a  wide  variety  of  problems  involving  both  complex  velocities  and  complex  p.  For  examples, 

discussion,  and  published  subroutine,  see  Cormier  and  Richards  (1989). 

* 

Although  individual  terms  in  a  ray  expansion  of  the  layered  response  function  have  poles  oriented 
at  ±60  degrees  with  respect  to  the  real  p  or  k  axis  (Scholte,  1956;  Nussenzweig,  1969;  Ludwig,  1970; 
Richards,  1973),  the  unexpanded  response  function  has  poles  only  on  the  real  p  or  k  axis.  A  simple 
illustration  of  how  the  ±60  degree  oriented  poles  vanish  in  the  complete  response  function  can  be 
made  by  considering  the  dispersion  function  for  Love  modes  in  a  vertically  inhomogeneous  layer 
overlying  a  high  velocity  cap  layer.  Substituting  SH  fundamental  matrices  from  AI.5  into  the 
propagator  matrix  in  AIL  15  gives  an  expression  for  the  Love  wave  dispersion  function: 


D1(Re)  =  +  *a?Ma)]ZMrc) 

+  *  PitoPifa  [»?1  rn  ejpe?  -  mm  frf^pafo)  (AII.17) 

where  the  subscripts  (1,2)  on  the  generalized  vertical  wavefunctions  and  slownesses  refer  to  evalua¬ 
tion  at  the  top  (r  =  Re)  and  bottom,  respectively,  of  the  layer  overlying  the  high  velocity  halfspace. 
The  Langer  approximation  is  a  uniform  asymptotic  approximation  of  fj,  i],  <r^\  cr^\  returning  the 
WKBJ  approximation  at  large  values  of  |wr|,  the  Airy  approximation  at  small  values  of  |«r|,  and 
a  smoQth  transition  between  the  WKBJ  and  Airy  approximation  at  intermediate  values  of  |wr|.  At 
large  values  of  jur|,  where  the  Langer  approximation  is  equivalent  to  the  WKBJ  approximation,  the 
±60  degree  oriented  poles  of  individual  reflection/ transmission /conversion  coefficients  are  oriented 
at  ±90  degrees  with  respect  to  tV  real  axis.  From  the  uniformity  of  the  Langer  approximation,  to 
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show  that  the  ±60  degree  poles  vanish  in  the  full  response  function,  it  is  sufficient  to  show  that  the 
±90  degree  poles  of  the  WKBJ  approximated  ray  expansion  vanish  in  the  full  response  function. 
Two  p  domains  of  behavior  should  be  considered,  corresponding  to  traveling  wavefunctions  and 
to  exponentially  decaying  and  growing  wavefunctions.  In  the  traveling  wave  domain,  the  WKBJ 
approximation  of  AII.17  becomes 


tan(n  - rj)  = (AIU8) 

This  is  identical  to  the  well  known  expression  for  the  dispersion  of  Love  waves  in  a  layer  over  a 
halfspace  if  (t\  -  r2)  is  replaced  by  the  vertical  slowness  in  the  homogeneous  layer  times  the  layer 
thickness  II,  i.e.,  (r\  -  r2)  — ►  \fp'^  -  ft  II.  Zeros  of  the  dispersion  function  AII.17  occur  at  real 
values  of  r  and  p. 

In  the  exponentially  decaying  p  domain,  no  traveling  waves  exist  in  the  layer  and  the  dispersion 
function  defined  by  AII.17  vanishes.  If  this  happens  during  the  upward  propagation  of  (Di,  112), 
the  calculation  of  D;  elements  can  be  started  at  a  higher  layer  (see  Layer  Reduction  subsection). 
This  behavior  in  the  exponentially  decaying  and  growing  domain  follows  from  the  properties  of  the 
generalized  vertical  wavefunctions  and  slownesses  at  large  values  of  \ut\  in  this  domain.  Specifically, 
for  large  values  of  |cvr|  in  the  exponentially  decaying  and  growing  domain  of  complex  p,  the  Love 
wave  dispersion  function  vanishes  because  the  Langer  subroutine  returns  77  =  -i?  and  c^1);  =; 
(Richards,  1973). 

In  the  p  domain  corresponding  to  rays  bottoming  below  the  surface  but  above  the  lower  bound- 
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ary,  the  dispersion  function  AII.17  may  be  rewritten  as: 


m  = 


P2% 


(AII.19) 


where  common  terms  have  been  canceled,  WKBJ  limits  of  the  vertical  wavefunctions  have  been 
taken  at  the  lower  boundary,  and  Langer  approximated  vertical  wavefunctions  at  the  surface  have 
been  recombined  using  the  property  that  -  r/<rW  =  The  generalized  vertical  slowness 

i)  is  defined  using  Airy  functions  of  the  type  Ai(-Q)  (see  AI.2  and  AI.l),  which  have  zeros  only 
on  the  positive  £p  axis.  Poles  of  the  dispersion  function  will  be  located  close  to  these  real  zeros. 
Since  (p  is  real  only  for  real  values  of  p,  the  poles  of  the  dispersion  function  must  lie  along  the  real 
p  axis. 


At  the  real  ray  parameters  corresponding  to  the  reciprocal  velocities  in  the  capping  layer,  branch 
cuts  are  oriented  along  the  real  p  axis,  making  a  right  angle  turn  at  the  p  origin  toward  the  positive 
imaginary  p  axis.  If  intrinsic  attenuation  is  assumed  with  complex  velocities  given  by  equations  4 
and  5,  the  poles  along  the  real  axis  are  shifted  slightly  into  the  upper  half  p  plane  when  the  inverse 
Fourier  time  transform  is  of  the  form  f(t)  =  /  F(u)e~,wt. 

The  integration  contours  leading  to  the  locked  mode  representation  are  nearly  identical  to  those 
described  and  shown  in  Harvey  (1981),  except  for  the  fact  that  the  contour  must  be  closed  in  the 
upper  rather  than  lower  half  of  the  p  plane.  The  difference  in  contour  closure  follows  from  a 


difference  in  the  sign  convention  chosen  for  the  Fourier  time  transform,  which  is  also  reflected  in  a 
difference  in  the  sign  of  the  imaginary  part  of  complex  velocity  and  the  use  of  7/W  instead  of  II W 
for  the  horizontal  wavefunction.  The  scalar  amplitude  functions  and  eigenfunctions  are  all  even 
functions  of  wavenumber  k  or  ray  parameter  p,  and  hence,  are  independent  cf  the  sign  convention 


used  for  Fourier  transforms.  This  makes  it  possible  to  use  Harvey’s  programs  for  the  construction 
of  complex  spectra  and  inversion  into  the  time  domain  simply  by  changing  the  sign  of  the  imaginary 
part  of  pole  positions  found  in  attenuative  structure. 


81 


Figure  1:  Discrete  (above)  and  continuous  (below)  representations  of  a  gradient  in  P  velocity  in  a 
test  model  of  the  crust. 


Figure  2:  Test  models  having  three  different  intensities  of  gradients  in  an  inhomogeneous  layer 
overlying  a  homogeneous  halfspace.  Model  1  is  the  test  model  of  Spudich  and  Ascher  (1983). 


Figure  3:  (a)  Love  and  (b)  Rayleigh  mode  dispersion  curves  calculated  in  Model  1  using  a  thin 
layered  representation  of  the  gradient  layer  (solid)  and  the  Langer  approximation  in  a  continuous 
representation  of  the  gradient  layer  (dashed). 


Figure  4:  Scalar  amplitude  factor  of  the  fundamental  Rayleigh  mode  plotted  against  the  nondi- 
mensional  parameter  of  wavelength/model  scale  length  in  the  models  1,  2,  and  3  shown  in  Figure 
2.  The  solid  curves  were  calculated  by  representing  the  surface  gradient  layer  with  40  thin  ho¬ 
mogeneous  layers;  the  dashed  curves  were  calculated  using  the  Langer  approximation  in  a  thick, 
continuous  surface  gradient  layer 


Figure  5:  Eigenfunction  for  vertical  displacement,  Eu  for  the  first  higher  Rayleigh  mode  at  0.5 
Hz  in  models  1,  2,  and  3  shown  in  Figure  2.  The  non-dimensional  parameter  X/L  was  calculated 
using  the  mean  shear  velocity  and  shear  velocity  gradient  in  the  surface  gradient  layer,  c  and  c’ 
are  respectively  the  phase  velocities  calculated  using  the  thin  homogeneous  layer  parameterization 
of  the  gradient  layer  and  a  thick  continuous  gradient  layer  with  the  Langer  approximation. 


Figure  6:  Comparison  of  synthetic  seismograms  calculated  in  Model  1  using  a  thin  layered  represen¬ 
tation  of  the  gradient  layer  (solid)  and  the  Langer  approximation  in  a  continuous  representation  of 
the  gradient  layer  (dashed).  The  source  is  a  point  double  couple  at  4.92  km.  depth,  corresponding 
to  a  vertically  dipping  strike  slip  fault,  striking  to  the  north,  observed  at  receivers  at  45°  azimuth.  A 
step  function  time  dependence  of  the  scalar  moment  is  assumed.  Shown  are  the  three  components 
of  particle  velocity.  The  effects  of  geometric  spreading  of  body  waves  have  been  approximately 
removed  by  multiplying  each  seismogram  by  range. 


Figure  7:  Comparison  of  synthetic  seismograms  calculated  in  Model  3  using  a  thin  layered  repre¬ 
sentation  of  the  gradient  layer  and  the  Langer  approximation  in  a  continuous  representation  of  the 
gradient  layer.  The  result  of  the  discrete  method  is  shown  at  each  range.  The  lower  amplitude 
trace  labeled  DIF  is  the  difference  between  the  seismograms  calculated  by  the  two  different  param- 
etcrizations,  (D)  discrete  thin  layered  and  (CL)  continuous  with  the  Langer  approximation,  i.c., 
DIF(t)  =  So(t)  -  Scl(1)‘  An  approximate  correction  for  geometric  spreading  of  body  waves 
has  been  made. 
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Figure  8:  Seismograms  for  the  fundamental  Rayleigh  mode  were  synthesized  in  a  test  anelastic 
mode),  (a)  Left:  surface  normalized  displacement  at  1  Hz  (solid)  and  0.1  Hz  (dashed),  (b)  Middle: 
P  and  S  velocity  at  1  Hz  (solid)  and  0.1  Hz  (dashed),  (c)  Right:  shear  attenuation,  Qjj1,  at  1  Hz 
(solid)  and  0.1  Hz  (dashed). 


Figure  9:  A  comparison  of  synthetics  for  the  vertical  component  of  the  fundamental  mode  Rayleigh 
wave  using  perturbation  theory  and  an  exact,  complex  mode  calculation.  At  each  range,  the  results 
of  the  exact  calculation  are  followed  by  the  differential  seismogram  obtained  by  subtracting  the 
seismogram  calculated  by  perturbation  theory  from  the  seismogram  calculated  by  complex  modes 
and  eigenfunctions.  Each  trace  is  normalized  by  its  peak  amplitude,  indicated  by  the  number  to 
the  left  of  each  trace. 


Figure  10:  A  simple  crust  and  upper  mantle  model  MH  composed  of  two  homogeneous  crustal  layers 
overlying  a  homogeneous  mantle  in  a  flat  earth.  Densities  and  attenuation  of  MH  are  those  given 
by  W-Y.  Kim  (1987).  (Note  slight  negative  gradients  in  MH  plotted  against  depth  in  a  spherical 
earth.)  Model  MG  has  a  single  crustal  gradient  layer  and  mantle  gradient  layer.  Densities  in  M1I 
assume  linear  gradients  in  depth  with  p  =  2.7,  2.9,  3.1,  and  3.8  at  depths  0,  38  km,  38  km,  and 
738  km.  Attenuation  in  MG  is  assumed  to  be  a  relaxation  band  in  pure  shear  between  0.0001  and 
4.0  Hz,  assuming  linear  gradients  with  depth  in  Q~l,  with  =  0.003,  0.0025,  0.005,  and  0.005 
at  0,  38  km,  38  km,  and  738  km.  Calculations  with  the  Langer-locked  mode  assumed  continuous 
velocity  functions  in  each  gradient  layer,  but  fine  scale  layering  was  used  to  simulate  the  gradients 
of  MG  in  the  wavenumber  integration  method. 


Figure  11:  A  comparison  of  synthetics  in  model  MH  (above)  and  MG  (below)  computed  by  the 
Langer-locked  mode  method,  summing  10  higher  Rayleigh  modes.  Shown  is  the  vertical  displace¬ 
ment  for  a  double  couple  point  source  at  30  km  depth.  The  orientation  of  the  double  couple 
corresponds  to  a  vertically  dipping  strike  slip  fault,  striking  to  the  north,  observed  at  an  azimuth 
of  45°.  A  step  function  ume  dependence  of  the  scalar  moment  is  assumed,  and  the  result  has  been 
convolved  with  a  short  period  WWSSN  instrument  response. 


Figure  12:  A  comparison  of  synthetics  in  model  MH  computed  by  the  Langer-locked  mode  method, 
summing  all  of  the  Rayleigh  modes  in  a  frequency  band  up  to  2  Hz  (above),  with  synthetics  in 
model  MH  computed  by  wavenumber  integration  (below). 


Figure  13:  A  comparison  of  synthetics  in  model  MG  computed  by  the  Langer  Langer-locked  mode 
method,  summing  all  of  the  Rayleigh  modes  in  a  frequency  band  up  to  2  Hz  (above),  with  synthetics 
in  model  MG  computed  by  wavenumber  integration  and  parameterization  of  gradient  layers  by  thin 
homogeneous  layers  (below). 
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